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DEVELOPMENT  OF  A  COMPREHENSIVE  MAGNETOHYDROOYNAMIC 
MODEL  OF  SOLAR-TERRESTRIAL  INTERACTION 


1.  EXECUTIVE  SUMMARY 

A  final  summary  report  is  provided  of  the  work  performed  under  AFOSR 
Contract  No.  F49620-86-C-0035.  The  work  under  this  contract  relates  to 
the  preliminary  development  of  a  comprehensive,  multi-level,  3-D 
magnetohydrodynamic  model  that  seeks  to  describe  the  detailed,  global 
interaction  of  the  solar  wind  plasma  and  magnetic  field  with  the  Earth’s 
geospace.  The  overall  technical  objective  of  the  present  work  is  the 
development  of  a  model  capable  of  quantitatively  describing  and 
routinely  predicting  many  of  the  most  important  space  plasma  and  field 
phenomena  associated  with  the  interaction  of  the  solar  wind  plasma  with 
the  terrestrial  environment.  The  ultimate  importance  of  such  a  model 
for  Air  Force  missions  in  the  Earth’s  geospace  lies  in  its  capability  to 
function  as  an  accurate  and  reliable  predictor  of  the  solar-terrestrial 
plasma  weather  and.  most  importantly,  its  ability  to  rapidly  warn  of 
impending  major  changes. 

Because  of  the  complexity  of  such  solutions  and  the  detail  required, 
the  methodology  for  such  a  model  must  necessarily  be  computational.  In 
the  current  work,  we  have  successfully  accomplished  the  development  of 
five  separate  advanced  computational  submodels  of  the  complete 
interaction  model  involving  approximately  42,000  lines  of  fortran  source 
code.  These  submodels  in  total  provide  the  capability  for 
quantitatively  simulating  a  variety  of  steady  and  unsteady  phenomena 
associated  with  the  solar  wind-terrestrial  environment  interaction 
process  to  a  degree  previously  unattainable.  A  number  of  these  new 
developments  have  been  incorporated  into  a  core  interaction  submodel 
that  can  be  structured  to  work  in  both  a  rapid  warning  mode  for 
spacecraft  protection  as  well  as  in  a  detailed  scientific  analysis  mode 
for  fundamental  studies.  The  other  submodels  developed  either  simulate 
similar  phenomena  of  the  solar  wind-terrestrial  environment  interaction 
process  as  the  core  interaction  submodel  but  at  a  higher  level  of 
simulation  accuracy  or  simulate  additional  interaction  phenomena. 

All  of  the  important  results  from  the  research  performed  under  this 
contract  have  been  reported  in  the  open  literature,  both  in  scientific 
archival  journals  and  as  technical  papers  at  scientific  meetings  with 
appropriate  acknowledgement  of  AFOSR  support.  This  report  provides  a 
synopsis  of  all  the  model  developments  accomplished,  a  summary  of  the 
important  research  results,  a  reference  list  of  the  publications  and 
presentations  resulting  from  the  research,  and  copies  of  all  the 
publ ications. 
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2.  INTRODUCTION 


The  space  era  has  provided  great  advances  in  the  understanding  of 
the  many  and  varied  phenomena  involved  in  solar  terrestrial  physics,  in 
particular,  those  involving  plasma  and  magnetic  fields  of  solar  origin. 
With  the  numerous  manned  and  unmanned  Air  Force  missions  currently  being 
planned  for  the  near-Earth  space  environment,  the  need  has  become  quite 
evident  for  procedures  to  describe  and  predict  plasma  and  field 
phenomena  associated  with  the  coupled  solar-interplanetary  medium- 
terrestrial  environment  system.  Essentially  all  of  these  future  space 
missions  will  require  a  reliable  on-going  forecast  of  solar-terrestrial 
plasma  "weather"  and  warning  of  impending  major  changes,  both  for  better 
utilization  of  observational  opportunities  as  well  as  for  the  safety  and 
protection  of  spacecraft.  In  addition,  many  important  terrestrial -based 
activities  such  as  communications  and  other  atmospheric  dependent 
applications  are  also  strongly  affected  by  high  levels  of  solar  plasma 
and  magnetic  field  activity.  As  the  space  station  and  associated  space 
activities  are  initiated  in  the  near  future,  the  need  for  such  a 
capability  will  become  even  more  critical. 

Although  there  have  been  numerous  and  significant  advances  in  our 
ability  to  predict  various  segments  of  the  coupled  solar-interplanetary 
medium-terrestrial  environment  system,  the  overall  quantitative  predic¬ 
tive  capability  presently  remains  limited  and  disjointed.  While  there 
currently  are  numerous  models,  some  highly  quantitative  and  verified  by 
comparisons  with  observations,  that  deal  with  various  components  of  the 
total  system,  these  models  normally  exist  as  separate  entities  and  have 
not  been  unified  with  others  to  any  significant  degree.  Additionally, 
even  fewer  models  deal  with  their  component  processes  with  the 
generality  required  to  make  the  prediction  usually  desired,  particularly 
with  respect  to  transient  phenomena. 

Additionally  and  perhaps  of  most  importance,  the  major  advancements 
that  have  been  achieved  in  fluid  dynamic  computational  methods  in  the 
past  15  years  have  not  been  applied  to  any  significant  degree  to  space 
physics  problems.  Because  most  of  the  important  problems  associated 
with  the  solar-terrestrial  environment  interaction  concern  plasma 
phenomena  that  are  complex,  interrelated,  and  time-dependent,  and 
because  the  solutions  for  these  phenomena  require  both  a  fine  detail  for 
resolution  as  well  as  a  global  extent,  the  ultimate  solution  methods 
must  necessarily  be  computationally  based.  A  primary  objective  of  this 
study  has  been  the  introduction  of  new  and  advanced  computational 
development  and  numerical  modeling  to  a  variety  of  detailed  and  advanced 
aspects  of  the  3-D  solar-terrestrial  interaction  system.  A  unique 
feature  of  our  model  is  that  certain  segments  contain  a  hierarchy  of 
computational  submodels  for  solving  the  same  plasma  physics  problem  but 
at  different  accuracy-of-simulation  levels.  This  multi-level  hierarchy 
enables  use  of  our  model  in  at  least  two  different  modes:  (1)  a  rapid 
analysis  mode  necessary  for  use  in  a  warning  or  preliminary  analyses 
capability,  and  (2)  a  detailed  analysis  mode  for  use  in  in-depth  scien¬ 
tific  studies.  This  feature  enables  accurate  simulation  of  a  variety  of 
complex  nonlinear  field  and  plasma  phenomena,  while  simultaneously 
providing  the  means  to  control  the  potentially  prohibitive  computer 
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costs  associated  with  detailed  studies  of  3-D  solar-terrestrial 
phenomena. 

During  this  phase  of  the  investigation,  we  have  successfully 
developed  several  completely  new  computational  submodels  for  simulating 
different  aspects  of  the  solar  wind  interaction  problem.  In  addition, 
we  have  also  developed  a  number  of  significant  extensions  and  enhance¬ 
ments  to  our  previously-developed  solar  wind/terrestrial  interaction 
model  (Refs.  1-3).  These  developments,  taken  in  total,  now  allow  us  to 
quantitatively  predict  details  of  the  plasma  and  field  properties 
associated  with  a  variety  of  interaction  phenomen  ,  both  steady  and 
unsteady,  to  a  degree  of  accuracy  and  routineness  that  was  heretofore 
unachievable.  These  modeling  developments,  which  are  described  in 
detail  in  the  following  sections,  have  been  Incorporated  into  a  core 
interaction  model  that  even  at  this  preliminary  stage  of  development 
establishes  a  comprehensive,  practical,  user-oriented,  working  model  for 
ourselves  and  other  interested  space  scientists  to  study  various  aspects 
of  the  global  3-D  solar  wind/terrestrial  interaction. 
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3.  TECHNICAL  OVERVIEW 


> 


> 


3.1  Overview  of  Solar-Terrestrial  System 

The  physical  foundation  of  the  ideas,  morphology  and  methodology 
underlying  the  current  development  program  is  illustrated  in  Figure  1 
below,  where  the  three  primary  segments  of  the  coupled  solar-terrestrial 
system  are  identified: 


THt  SOLAR  -  TERRESTRIAL  SYSTEM 


I 


Solar  1  Interplanetary 
Initiation  i  Propagation 

Segment  I  Segment 


Figure  1.  Illustration  of  the  Three  Primary  Segments  of  the  Coupled 
Solar-Terrestrial  System 


The  outer  layers  of  the  sun  are  a  convection  heat  engine  that 
produce  both  large-  and  small-scale  hydromagnetic  motions  on  the  solar 
surface.  This  activity  at  the  solar  surface  creates  the  first  segment 
of  the  solar-terrestrial  system  by  generating  the  solar  wind--the 
underlying  medium  that  couples  the  solar-terrestial  system.  The  solar 
wind  itself  arises  from  the  hydromagnetic  motions  on  the  solar  surface 
which  provide  a  complex  magnetic  topology  from  which  the  heated  solar 
plasma  ultimately  escapes  into  interplanetary  space.  The  solar  wind 
plasma  thus  created  carries  a  portion  of  the  solar  magnetic  field  as  it 
subsequently  flows  from  the  solar  corona  out  into  the  solar  system.  The 
solar  wind  itself  is  a  supersonic,  super-Alfvenic.  strongly  ionized 
plasma  that  transports  plasma,  energy,  angular  momentum,  and  magnetic 
field  from  the  sun  past  all  the  planets  of  the  solar  system.  The  second 
segment  of  the  solar-terrestrial  system  is  the  interplanetary 
propagation  segment.  It  involves  the  transport  through  the  solar  wind 
media  of  both  steady  and  quasi-steady  processes  commonly  occurring  in 
the  solar  wind,  as  well  as  the  propagation  of  high  energy  transient 
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phenomena  such  as  solar  flare  Initiated  disturbances.  This  segment 
initiates  at  the  solar  corona  and  and.  for  our  purposes,  ends  in  the 
near-vicinity  of  the  Earth's  geospace.  The  final  segment  of  the  coupled 
solar-terrestrial  system  is  the  terrestrial  response  segment.  This 
segment  involves  the  global  interaction  and  subsequent  response  of  the 
Earth's  environment  to  both  the  steady  state  and  transient  plasma  and 
field  phenomena  that  have  been  transported  and  modified  by  the  solar 
wind  plasma  from  the  surface  of  the  sun.  through  the  interplanetary 
medium,  and  finally  impact  on  the  Earth's  geospace. 

Each  of  the  three  segments  encompasses  varied  and  complex  plasma  and 
magnetic  field  phenomena,  and  each  in  itself  actually  represents  an 
important  subdiscipline  of  solar-system  plasma  physics.  The  ultimate 
predictive  capability  for  the  coupled  solar-terrestrial  system  will 
involve  detailed  predictive  models  for  all  three  segments  integrated 
with  each  other  and  coupled  with  simultaneous  observational  input  data 
provided  by  multiple  spacecraft  in  the  environments  of  the  near-sun. 
interplanetary  space,  and  near-Earth.  The  focus  for  the  current  study 
is  on  the  terrestrial  response  segment  and  represents  the  logical  first 
step  in  achieving  the  complete  predictive  capability  for  the  sun- 
terrestrial  system.  The  initiating  input  data  for  the  interaction 
models  developed  here  is  assumed  to  be  provided  by  one  or  more 
monitoring  spacecraft  located  some  distance  upstream  of  the  Earth  toward 
the  sun.  That  observational  data  is  then  projected  via  use  of  the 
interaction  model  to  predict  plasma  and  field  properties  from  the 
spacecraft  monitor  location  through  the  bow  shock  and  magnetosheath 
region  down  to  the  terrestrial  magnetopause  boundary. 


3.2  Global  Interaction  Submodels  of  Total  Interaction  Model 

Because  of  the  many  and  varied  plasma  and  field  phenomena  that  we 
intend  to  describe,  it  is  neither  convenient  nor  efficient  for  the 
terrestrial  interaction  model  to  be  composed  of  a  single,  large,  all- 
encompassing  computational  element.  We  have  chosen  instead  to  configure 
the  total  model  as  a  modularized  assemblage  of  individual  submodels  that 
are  linked  together  to  form  the  complete  interaction  model.  In  addition 
to  this  modularization,  we  have  also  intentionally  constructed  the  model 
to  contain  a  hierarchy  of  different  accuracy-of-simul ation  levels  and 
associated  computational  algorithms  in  order  to  simulate,  when  necessary 
and  at  the  appropriate  level  of  simulation  accuracy,  specific  plasma  and 
magnetic  field  phenomena  associated  with  the  global  solar  wind- 
terrestrial  magneosphere  interaction  process. 

In  broad  summary,  the  investigation  completed  here  has  resulted  in  a 
comprehensive  interaction  model  in  preliminary  form  that  provides,  on  a 
single-fluid  continuum  rather  than  multi -component  particle  basis,  a 
quantitative  three-dimensional  description  of  the  spatially  and 
temporally  varying  properties  of  the  plasma  and  magnetic  fields  that 
link  the  near-Earth  interplanetary  medium  with  the  terrestrial 
magnetopause.  The  complete  interaction  model  is  capable  of  predicting 
the  global  plasma  and  magnetic  field  properties  associated  with  the 
transport  of  the  interplanetary  solar  wind  plasma  from  an  arbitrary 
location  of  a  spacecraft  monitor  placed  upstream  of  the  terrestrial  bow 
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shock  to  and  through  the  bow  shock  and  intervening  magnetosheath  region 
and  down  to  the  magnetopause  boundary.  The  model  consists  of  several 
distinct  but  interactive  modules  that  quantitatively  predict  various 
steady  and  dynamic  aspects  of  the  interaction.  The  model  employs 
advanced  computational  procedures,  many  of  which  are  of  recent 
development. 

In  a  conceptual  but  directly  realizable  operational  mode,  the  entire 
model  calculation  can  be  configured  to  be  carried  out  at  a  pace  that  is 
at  least  two  orders  of  magnitude  faster  than  the  actual  physical  events 
evolve  in  space.  We  view  such  an  operational  mode  of  the  final 
interaction  model  as  a  pilot  system  capable  of  providing  a  solar- 
terrestrial  plasma  weather  prediction  service  analogous  partly  to  the 
global  weather  forecasting  system,  and  partly  to  the  Pacific  Ocean 
tsunami  warning  system.  The  complete  implementation  of  the  prediction 
and  warning  mode  of  the  model  will  require  that  the  interplanetary 
inputs  be  continually  monitored  and  provided  in  an  on-line  mode  to  the 
interaction  model.  Both  the  slowly  varying  pseudo-steady  state  of  the 
interplanetary  medium  and  the  propagation  of  disturbances  through  it  can 
then  be  simulated  and  routinely  and  accurately  predicted  via  use  of  the 
interaction  model . 

In  order  to  realize  such  a  quantitative  predictive  capability  of  the 
global  interaction  process,  we  have  proceeded  in  this  contract  to  both 
extend  previously  existing  interaction  submodels  and  to  develop  entirely 
new  interaction  submodels.  The  effort  has  involved  the  successful 
development  and  implementation  of  several  important  extensions  in  the 
submodel  that  we  have  termed  the  basic  or  core  interaction  model,  in 
particular  the  prediction  of  time-dependent  dynamic  processes  involving 
the  interplanetary  magnetic  field.  In  additional,  we  have  successfully 
developed  a  number  of  new  interaction  submodels.  During  this  contract 
effort,  four  completely  new  computational  submodels  have  been  developed 
which  provide  predictive  capabilities  beyond  the  basic  core  interaction 
model.  Each  of  these  new  submodels  employs  different  advanced  numerical 
algorithms,  and  solves  in  detail  a  different  aspect  of  the  solar  wind- 
terrestrial  environment  interaction  problem.  These  new  interaction 
submodels  are: 


•  3-D  steady  flow  interaction  model 

•  Aligned  flow  full  MHD  interaction  model 

•  Continuous  flow  transient  interaction  model 

•  Discontinuous  flow  transient  interaction  model 


The  developments  associated  with  these  new  submodels  as  well  as  those  of 
the  core  interaction  model  are  discussed  in  detail  in  the  following 
section. 
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4. 


SUMMARY  OF  INTERACTION  MOOEL  DEVELOPMENTS 


4.1  Development  of  Core  Interaction  Model 

Central  to  the  predictive  model  being  developed  here  is  the 
submodel  that  we  have  termed  the  core  interaction  model.  The  core 
interaction  model  represents  the  basic  or  first-order  computational 
model  of  the  complete  multi-level  interaction  model  for  predicting  the 
global  3-D  interaction  of  the  solar  wind  with  the  Earth's  magnetosphere. 
For  the  current  program,  the  core  model  provides  the  fundamental 
computational  tool  for  determining  the  detailed  plasma  and  field 
properties  in  the  magnetosheath  region  between  the  Earth's  bow  shock  and 
magnetopause. 

The  theoretical  and  computational  basis  of  the  current  model  is 
described  in  detail  in  Refs.  1-3.  Here  we  present  only  some  salient 
features  of  the  model  convenient  for  the  present  discussion.  The  funda¬ 
mental  assumption  of  the  model  is  that  the  bulk  properties  of  solar  wind 
flow  past  planetary  magnetoi onopause  obstacles  can  be  adequately 
described  by  solutions  of  the  continuum  equations  of  magneto¬ 
hydrodynamics  for  a  single-component  perfect  gas  having  infinite  elec¬ 
trical  conductivity  and  zero  viscosity  and  thermal  conductivity.  The 
governing  equations  express  the  conservation  of  mass,  momentum,  energy, 
and  magnetic  field,  augmented  by  jump  conditions  appropriate  for  a  fast 
shock  wave  and  a  tangential  discontinuity  at.  respectively,  the  bow 
shock  wave  and  magnetoi onopause  surface.  The  model  determines  the 
magnetoionopause  and  bow  shock  locations  and  the  magnetosheath  plasma 
and  field  properties.  Justification  of  the  theoretical  basis  rests 
primarily  on  the  outstanding  agreement  between  results  calculated  in 
this  way  and  observations  made  in  space  for  the  Earth  and  other  planets 
over  the  past  25  years.  This  basis  is  now  well  established. 

The  basic  or  core  interaction  model  does  not  solve  the  general 
unsteady  3-D  MHD  equations  but  a  simplified  subset  based  on: 


•  steady  flow 

•  predetermined  magnetopause  obstacle 

•  high  Alfvenic  Mach  number 


The  last  assumption,  typical  of  solar  wind  flows  past  the  planets, 
results  in  the  magnetic  terms  in  the  momentum  and  energy  equations  being 
significantly  smaller  than  the  gas  dynamic  terms.  This  implies  that  the 
solution  for  the  fluid  motion  can  be  decoupled  from  the  magnetic  field. 
The  simplified  equations  governing  the  flow  field  then  reduce  to  the 
familiar  gas  dynamic  Euler  equations.  The  magnetic  field  can  then  be 
determined  subsequently  based  upon  the  calculated  flow  field.  It  has 
been  shown  that  the  magnetic  field  determined  in  this  manner--when  using 
the  exact  MHD  rather  than  approximate  gas  dynamic  flow  field--obeys  the 
conservation  principle  that  the  magnetic  flux  passing  through  an  arbi¬ 
trary  surface  moving  with  the  fluid  is  conserved  or,  in  more  picturesque 
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terms,  ;.hat  the  magnetic  field  is  convected  with  the  fluid.  It  is 
generally  regarded  that  the  magnetic  field  calculated  in  this  way  i.e.. 
when  employing  the  approximate  gas  dynamic  flow  field,  accurately 
predicts  the  general  features  as  well  as  many  details  of  the  observa¬ 
tions.  This  has  been  overwhelmingly  confirmed  with  the  current  model  by 
the  many  comparisons  now  made  of  convected  field  results  with 
observations  at  Earth.  Venus,  and  Mars  (see  Refs.  4-8).  The 
Interaction  solution  determined  in  this  fashion  is  now  commonly  referred 
to  as  the  gas  dynamic-convected  field  approximation. 

The  decoupling  of  the  B  field  calculation  from  the  flow  field 
determination  is  a  particularly  convenient  and  valuable  feature  in  the 
basic  interaction  model.  The  flow  field  determination  is  the  most  time 
consuming  portion  (-90%  CPU  time)  of  the  total  computational  process  for 
a  single  case.  However,  once  the  flow  field  is  determined  and  stored, 
many  different  magnetosheath  fields  due,  for  example,  to  different  IMF's 
can  be  quickly  and  inexpensively  determined  as  a  separate  process. 
Because  the  B  calculation  is  extremely  rapid,  this  feature,  now  imbedded 
in  the  model,  has  been  exploited  in  a  number  of  recent  investigations, 
for  example  Refs.  4-8. 

In  order  to  determine  the  flow  field,  two  separate  gas  dynamic 
solvers,  as  shown  in  Figure  2  below,  have  been  incorporated  in  the  model 
and  linked  together. 
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Nose  Region  |  Tail  Region 


Unsteady  gasdynamic  Steady  gasdynamic 
time-marching  solver  spatial-marching 


Axisymmetric  magnetoionospheric  obstacle 


Nose  Region  Solver 

•  Beam  and  Warming  implicit  algorithm 

•  Transient  and  asymptotic  steady  state  .olution  capabilities 

•  Fitted  bow  shock  and  obstacle  discontinuity  surfaces 

•  Fully  implicit  boundary  conditions 


Tail  Region  Solver 

•  Shiff  and  Steger  implicit  algorithm 

•  One  step  steady  state  solver 

•  Fitted  bow  shock  and  obstacle  discontinuity  surfaces 

•  Fully  implicit  boundary  conditions 


Figure  2.  Illustration  of  NOSE/TAIL  Subdivision  of  Flow  Field  and 
Associated  Solvers  Incorporated  into  Present  Core  Interaction  Model 


These  flow  solvers  are: 


NOSE  Region  Solver  -  to  determine  the  flow  field  in  the 
subsolar  region  and  sufficiently  far  downstream  to  encompass 
the  entire  subsonic  flow  portion  of  the  flow  field 


TAIL  Region  Solver  -  to  determine  the  flow  field  downstream  of 
the  nose  region  to  any  arbitrary  downstream  distance 
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The  core  model  has  the  capability  for  use  in  both  a  rapid  warning 
mode  to  warn  of  impending  changes  in  the  magnetosheath  plasma  and  field 
due  to  observed  solar  disturbances  propagating  in  the  interplanetary 
media,  and  also  for  use  in  a  more  detailed  scientific  analysis  mode  for 
fundamental  research  studies.  In  the  current  development  of  the  core 
model,  we  have  now  extended  and  greatly  enhanced  the  original  model  that 
was  previously  developed  (Refs.  1-3).  The  model  was  completely 
rewritten  to  incorporate  refined  computational  flow  and  magnetic  field 
solvers.  It  has  now  been  developed  into  a  completely  user-oriented 
structure,  and  has  been  extensively  tested  to  insure  robustness, 
generality,  and  ease  of  operation.  An  extensive  graphics  capability  has 
been  added  so  as  to  allow  the  complete  and  detailed  examination  of  all 
plasma  and  field  quantities  in  the  interaction  region.  Importantly,  the 
core  model  also  contains  the  capability  for  providing  time  histories  of 
the  model -predicted  plasma  and  field  properties  along  any  arbitrary 
spacecraft  trajectory  passing  through  the  magnetosheath  flow  field.  We 
and  others  have  found  this  particular  feature  to  be  invaluable  for 
comparative  studies  with  observational  data. 


Included  in  the  output  capability  of  the  core  interaction  model  are 
the  following: 


•  Contours  of: 

*  Plasma  temperature 

*  Plasma  pressure 

*  Plasma  density 

*  Plasma  velocity  magnitude 

*  Plasma  sonic  Mach  number 

*  Magnetic  field  components 

•  Plasma  flow  streamlines 

•  Plasma  flow  characteristic  lines 

•  Vector  plots  of: 

*  Plasma  velocity 

*  Magnetic  field 

•  Spacecraft  trajectory  plots: 

*  Ecliptic  plane  view 

*  Cylindrical  plane  view 

•  Time  histories  along  spacecraft  trajectory 

*  Plasma  pressure,  temperature,  density 

*  Velocity  magnitude  and  components 

*  Magnetic  field  magnitude  and  components 
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•  Diagnostic  analysis  plots  along  any  segment  of  the 
computational  grid: 

*  Plasma  pressure,  temperature,  density 

*  Velocity  magnitude  and  components 

*  Magnetic  field  magnitude  and  components 

A  brief  subset  of  the  detailed  graphics  output  is  provided  in 
Figures  3  and  4. 


Figure  3.  Sample  of  Selected  Graphic  Plasma  and  Field  Results 
Provided  by  the  Core  Interaction  Model 
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Figure  4.  Example  of  Plasma  and  Field  Results  Standardly  Provided  by 
the  Core  Model  as  Time  Histories  and  Vector  Plots  Along  a  Spacecraft 
Trajectory 


Figures  5  and  6  provide  examples  of  the  actual  use  of  such  results 
from  the  core  model  that  were  directly  employed  in  two  recent  space 
physics  applications  (Refs.  4,  5).  Figure  5  displays  a  comparison  of 
model -predi cted  magnetic  field  vectors  with  observations  along  a 
spacecraft  trajectory,  while  Figure  6  displays  a  prediction  of  locations 
of  magnetic  merging  sites  on  the  Earth’s  magnetopause  as  a  function  of 
IMF  orientation.  As  a  final  note,  the  source  code  for  the  graphics 
drivers  was  intentionally  written  by  us  rather  than  developed  from 
proprietary  routines.  This  was  done  to  enable  the  entire  package  of 
numerical  plus  graphics  routines  embodied  in  the  core  model  to  be 
directly  transportable  without  restriction  to  interested  Air  Force  space 
scientists  and  others. 


ORBIT  432 


Compared  With  Actual 


ORBIT  4  SB 
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el  Results  for  Magnetic  Field  Directly 
Spacecraft  Observations  (Ref.  4) 


Contours  of  squat  values  ol 
the  cosine  of  the  angle 
between  the  magnetospherfc 
and  the  magnetosheath 
magnetic  fields  for  various 
values  (B*.  By.Bz)  of  the  IMF 


Figure  6.  Example  of  Core  Model  Results  for  Magnetic  Merging  Site 
Locations  on  the  Earth's  Magnetopause  as  Viewed  From  the  Sun  for 
Various  IMF  Orientations  (Ref.  5) 


In  addition,  a  new  computational  module  has  now  been  developed  to 
extend  the  capability  of  the  core  model  to  predict  the  time-dependent 
magnetic  field  anywhere  in  the  global  interaction  region  due  to  an 
unsteady  IMF  propagating  past  the  magnetopause.  This  procedure  is 
capable  of  handling  both  a  continuously  varying  or  discontinuous  IMF. 
The  computational  submodel  employs  the  basic  core  model  for  the  flow 
field  and  applies  a  general  convected-field  calculation  procedure  for 
determining  the  unsteady  B  field  anywhere  in  the  magnetosheath.  Figure 
7  provides  an  example  of  this  capability  and  displays  the  predicted 
unsteady  magnetic  field  at  several  points  in  the  magnetosheath  due  to 
the  sinusoidally  time-varying  IMF  indicated.  In  addition  to  this  direct 
prediction  capability,  whereby  based  upon  a  given  unsteady  IMF  the 
resultant  unsteady  B  field  can  be  predicted  anywhere  in  the  magne¬ 
tosheath.  we  have  also  developed  a  computational  module  for  the  inverse 
problem.  In  this  mode,  if  the  unsteady  B  field  is  measured  by  a  space¬ 
craft  at  any  point  in  the  magnetosheath,  the  inverse  procedure  can 
determine  the  IMF  time  history.  This  provides  a  convenient  means  to 
estimate  the  unsteadiness  in  the  IMF  based  upon  single  spacecraft 
measurements,  and  a  direct  comparative  capability  of  prediction  and 
observation  for  a  two  spacecraft  system  where  one  spacecraft  serves  as 
an  upstream  monitor.  As  in  all  of  the  convected-field  computations  in 
the  core  model,  this  procedure  is  extremely  rapid  requiring  less  than  a 
minute  on  any  superminicomputer  system  for  an  entire  trajectory  computa¬ 
tion.  The  complete  core  interaction  model,  including  the  unsteady 
magnetic  field  prediction  module  and  the  display  graphics,  comprises 
approximately  20.000  lines  of  fortran  source  code. 
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Figure  7.  Illustration  of  the  Capability  of  the  Extended  Core 
Interaction  Model  to  Predict  Unsteady  Magnetic  Fields  in  the 
Magnetosheath  as  a  Function  of  a  Time-Varying  IMF 


4.2  Development  of  Full  3-D  Interaction  Model 

One  of  the  major  developments  of  the  modeling  effort  planned  for 
this  program  was  the  extension  of  the  numerical  solvers  embodied  in  the 
core  model,  which  apply  to  an  axisymmetric  magnetopause  obstacle 
geometry,  to  a  full  arbitrary  3-D  magnetopause  geometry  capability.  For 
reasons  of  computational  efficiency  as  pointed  out  in  the  previous 
section  discussing  the  core  interaction  model,  all  of  the  various 
interaction  models  developed  in  this  investigation  have  been  developed 
to  contain  two  separate  and  distinct  flow  field  solvers.  This 
particular  computational  feature  is  unique  to  the  current  model  in 
comparison  to  all  existing  MHD  models,  and  provides  a  number  of 
important  capabilities  in  our  model  that  are  unattainable  with  any  of 
the  other  existing  MHD  models.  The  first  flow  field  solver  is  an 
asymptotic  time-marching  solver  and  is  employed  to  determine  the  flow 
field  in  the  subsol ar-to-termi nator  region  (NOSE  solver),  while  the 
other  solver  spatially  marches  that  solution  downstream  from  the 
terminator  plane  as  far  as  required  (TAIL  solver). 

We  have  now  completed  the  development  of  both  of  the  flow  solvers 
for  the  3-D  model.  In  order  to  test  and  verify  these  solvers,  we  have 
performed  extensive  numerical  tests  of  the  3-D  model  on  several  classes 
of  3-D  magnetospheric  shapes.  Figure  8  provides  results  of  one  such 
application  that  was  related  to  a  study  of  the  effect  of  asymmetry  in 
the  Earth's  magnetopause  on  the  bow  shock  location  and  shape.  The 
insert  in  the  bottom  of  that  figure  provides  a  prospective  view  of  one 
quadrant  of  the  typical  computational  grid  employed  in  the  3-D  model. 
Illustrated  in  that  insert  are  the  complete  computational  grids  in  the 
equatorial  (X.Y)  and  polar  (X.Z)  azimuthal  planes,  together  with  the 
intersection  of  the  grids  in  the  other  azimuthal  planes  with  the  magne¬ 
topause  surface.  Shown  next  to  the  computational  grid  is  a  sketch 
illustrating  the  distortion  from  circular  symmetry  in  the  cross 
sectional  shape  of  the  Earth's  magnetopause  that  was  systematically 
varied  in  this  study.  The  three  plots  in  the  upper  part  of  Figure  8 
provide  the  calculated  3-D  model  results  for  the  bow  shock  locations  and 
the  corresponding  magnetopause  locations  in  the  10  equally-spaced 
azimuthal  planes  shown  in  the  computational  grid.  The  results  in  those 
plots  are  for  Mach  -  8.0  flow  past  a  terrestrial  magnetopause  shape 
whose  cross-section  has  been  distorted  from  a  circular  shape  into  an 
elliptic  shape  In  which  the  ellipse  has  its  major  axis  in  the  polar 
meridian  and  its  minor  axis  in  the  equatorial  plane.  The  plot  on  the 
upper  left  of  Figure  8  is  for  a  cross  sectional  shape  with  a  major/minor 
axis  ratio  A/B  -  1.05.  Corresponding  results  for  modeled  magnetopauses 
that  have  elliptic  cross  sections  with  major/minor  axis  ratios  of  A/B  - 
1.10  and  1.15  are  Illustrated  In  the  middle  and  right-hand  plots, 
respectively  for  these  various  elliptic  shapes.  The  two  plots  on  the 
bottom  left  and  right  of  Figure  8  display  a  summary  of  the  bow  shock  and 
magnetopause  locations  in  the  equatorial  and  polar  meridian  planes, 
respectively.  These  results  taken  In  total  show  that  the  typical  bow 
shock  surface  displacement  from  axisymmetry  that  occurs  due  to  asymmetry 
effects  in  the  magnetopause  surface  is  somewhat  less  than  approximately 
50*  of  the  corresponding  displacement  from  axisymmetry  of  the 
magnetopause,  an  important  result  to  know  when  trying  to  evaluate  or 
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ascertain  3-0  geometry  distortions  in  the  Earth's  bow  shock  and 
magnetopause  shapes. 


Figure  8.  Results  from  the  3-D  Terrestrial  Interaction  Model 
Illustrating  the  Effect  of  Asymmetry  in  the  Earth's  Magnetopause 
Cross  Sectional  Shape  on  the  Terrestrial  Bow  Shock  Shape  and 
Location 


Since  the  3-D  geometric  distortion  of  the  Earth’s  magnetopause  is 
not  as  severe  as  those  of  some  other  planets  in  the  solar  system,  we 
have  used  joint  NASA  support  to  provide  the  means  to  carry  out  a  more 
severe  test  and  verification  of  the  3-D  model.  Employing  observational 
data  obtained  for  Jupiter  and  Saturn,  we  have  used  the  3-D  model  to 
determine  the  flow  fields  and  corresponding  bow  shock  surfaces  that 
exist  about  these  highly  three-dimensional  obstacles.  The  computational 
model  performed  exceptionally  well  in  these  tests,  and  we  observed  no 
convergence  difficulties  for  all  of  the  extensive  numerical  testing 
performed.  Furthermore,  based  on  the  parametric  studies  that  were 
carried  out  involving  the  systematic  variation  of  the  3-D  magnetopause 
geometry,  the  3-D  predictive  model  was  able  to  quantitatively  determine 
the  amount  of  3-D  distortion  that  must  exist  in  these  magnetopauses  and. 
as  a  consequence,  resolve  the  anamoly  previously  noted  (Ref.  9)  with 
regard  to  the  observed  equatorial'  bow  shock  locations  for  these  planets. 


The  results  of  this  study  were  reported  in  Ref.  10.  and  were  also 
presented  at  the  Spring  1987  AGU  and  Summer  AIAA  meetings.  A  complete 
description  of  the  details  of  the  numerical  algorithms  employed  in  the 
3-D  model  is  provided  in  Ref.  10.  Finally,  a  revised  version  of  Ref.  10 
is  presently  in  review  for  publication  in  the  JGR. 
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4.3  Development  of  Aligned  Flow  Full  MHD  Interaction  Model 

The  final  steady  state  computational  interaction  submodel  that  has 
been  developed  extends  the  level  of  simulation  accuracy  beyond  the  gas 
dynamic-convected  field  approximation  employed  in  the  other  interaction 
submodels  to  that  of  full  MHD.  This  submodel  is  capable  of  predicting 
the  global  solar  wind-planetary  interaction  at  the  full  MHD  accuracy 
level  for  situations  when  the  IMF  is  aligned  with  the  oncoming  flow 
direction.  The  submodel  is  appropriate  for  determining  supersonic. 
super-Alfvenic  solar  wind  flows  past  general  axisymmetric  magnetosphere 
obstacle  shapes.  As  with  all  the  other  submodels  developed,  this 
computational  model  also  employs  two  separate  but  coupled  flow  field 
solvers  to  determine  the  steady  flow  field:  an  asymptotic  time-marching 
procedure  for  the  subsolar  to  terminator  region,  and  a  spatially- 
marching  procedure  for  post-terminator  locations.  Both  of  these  new 
computational  procedures  are  based  on  current,  fully-impl icit  numerical 
algorithms,  and  both  employ  fitted  discontinuity  representations  for  the 
bow  shock  and  magnetoi onospher i c  obstacle  surfaces.  Initial 
applications  of  the  model  have  now  been  successfully  made  to  Earth  ar.d 
Venus.  Some  results  of  application  of  the  new  aligned  flow  MHD  model  to 
Earth  are  provided  in  Figure  9.  In  that  figure,  results  from  the  new 
model  are  compared  to  a  selected  number  of  comparable  results  previously 
determined  by  Spreiter  and  Rizzi  (Ref.  11)  using  a  completely  different 
computational  method.  The  new  model  predicts  the  identical  variation  of 
the  bow  shock  shape  as  a  function  of  Alfven  Mach  number  in  the  oncoming 
solar  wind  that  was  previously  obtained,  that  is,  as  the  oncoming  Alfven 
Mach  number  decreases,  the  bow  shock  simultaneously  moves  inward  toward 
the  magnetopause  in  the  nose  region  and  flares  outward  away  from  the 
magnetopause  along  its  flanks.  Moreover,  it  attains  in  a  much  more 
effective  and  general  way  than  before  the  capability  for  routinely 
calculating  exact  MHD  solutions  for  aligned  flows  past  a  planetary 
magnetosphere.  Such  results  from  the  new  model  demonstrate  a  new 
capability  for  quantitatively  determining  the  sensitivity  of  the 
terrestrial  bow  shock  to  MHD  effects.  In  particular,  the  new  model 
provides  a  direct  means  for  evaluating  the  dependence  of  the  bow  shock 
location  and  shape  on  the  Alfven  Mach  number  of  the  oncoming 
interplanetary  solar  wind.  More  recent  studies  with  the  new  model  have 
unequivocally  demonstrated  the  parametric  importance  of  the  magnetosonic 
Mach  number  for  characterizing  MHD  effects  in  planetary  bow  shocks. 
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Figure  9.  Bow  Shock  Locations  Predicted  by  the  New  Aligned  Flow 
Submodel  Compared  With  Spreiter  and  Rizzi  (Ref.  11)  for  Aligned  Flow 
Past  the  Earth’s  Magnetopause  with  Sonic  Mach  Number  M«  -  10.  y  -  5/3, 
and  Various  Alfven  Mach  Numbers  Ma« 


4.4  Development  of  Continuous  Flow  Transient  Interaction  Model 

The  continuous  flow  transient  interaction  model  provides  the 
capability  for  accounting  for  transients  that  may  occur  in  the  oncoming 
interplanetary  plasma  flow  that  propagate  into  an  established  bow  shock 
and  magnetosheath  flow  field.  The  model  employs  two.  new,  fully- 
implicit  computational  solvers  coupled  in  the  familiar  NOSE/TAIL 
combination  employed  in  the  core  model.  Additionally,  a  new  computa¬ 
tional  module  has  been  added  to  determine  the  convected  B  field  in  a 
highly  refined  manner  employing  a  dense  streamline  coverage  and  a  grid 
clustering  capability.  This  important  feature  will  be  integrated  into 
the  core  interaction  model  in  the  near  future.  Figure  10  illustrates 
the  extended  upstream  grid  employed  in  the  new  model  to  account  for  the 
propagation  of  these  transients  into  the  bow  shock  and  magnetosheath, 
and  the  high  resolution  and  detail  in  the  predicted  flow  magnetic  field 
in  the  magnetosheath  region  that  is  achievable  from  this  model. 
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Figure  10.  Illustration  of  the  Extended  Upstream  Grid  and  High- 
Resolution  Plasma  and  Field  Predictive  Capability  of  the 
Continuous  Flow  Transient  Interaction  Submodel 
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The  significance  of  the  continuous  flow  transient  interaction  model 
to  this  development  effort  is  threefold:  first,  the  model  provides  an 
important  upgrade  of  the  computational  flow  models  employed  in  the  basic 
core  interaction  model;  second,  it  provides  the  capability  for  an 
extremely  refined  computation  of  the  plasma  and  convected  magnetic  field 
properties  in  the  vicinity  of  the  magnetopause  obstacle,  particularly  in 
the  subsolar  region;  and  third,  it  provides  an  accurate  numerical  model 
for  quantitatively  testing  simpler  rational  unsteady  interaction  models 
for  predicting  time-dependent  dynamic  plasma  properties  due  to  unsteady 
behavior  in  the  oncoming  solar  wind  plasma.  One  of  the  most  Important 
results  of  this  modeling  effort  has  been  the  unambiguous  demonstration 
that  simulations  of  unsteady  dynamic  processes  involving  the  interaction 
of  the  solar  wind  with  the  Earth’s  magnetopause  involve  so  much 
computational  resources  with  regard  to  computational  time  and  memory 
that  for  such  dynamic  interaction  models  to  be  useful  they  must  be 
computationally  simplified.  As  discussed  previously  with  regard  to  the 
core  interaction  model,  development  of  a  time-dependent  dynamic  model  to 
predict  the  behavior  of  the  convected  magnetic  field  based  upon  the 
underlying  plasma  flow  properties  remaining  steady  is  relatively 
straightforward,  and  has  now  been  accomplished  within  the  core 
interaction  model.  However,  the  corresponding  development  cf  a  time- 
dependent  dynamic  model  for  predicting  the  plasma  properties  in  the 
interaction  region  is  more  involved.  Nevertheless,  the  continuous  flow 
transient  interaction  model,  and  the  discontinuous  flow  transient 
interaction  model  discussed  in  the  following  section,  provide  the  means 
to  develop  and  quantitatively  verify  such  dynamic  models.  We  view  the 
development  of  more  efficient  time-den,-  Jent  plasma  and  magnetic  field 
models  as  the  next  most  importan*  step  in  modeling  the  solar  wind- 
terrestrial  interaction  process. 
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4.5  Development  of  Discontinuous  Flow  Transient  Interaction  Model 

The  discontinuous  flow  transient  interaction  model  is  capable  of 
predicting  the  unsteady  global  interaction  that  occurs  when  discontin¬ 
uous  transients,  such  as  interplanetary  shocks,  impinge  on  the  terres¬ 
trial  environment  and  subsequently  interact  with  an  established  bow 
shock  and  magnetopause  configuration.  This  interaction  submodel  also 
employs  an  extended  upstream  grid,  but  instead  of  using  a  fitted  discon¬ 
tinuity  representation  for  the  bow  wave,  as  is  used  in  both  the  core  and 
continuous  flow  transient  interaction  models,  this  model  uses  a  shock 
capturing  procedure  to  capture  all  shocks  in  the  flow  field.  This  is 
necessary  as  shock-on-shock  interactions  are  too  cumbersome  to  treat 
with  a  fitted  discontinuity  representation.  While  the  theoretical  basis 
underlying  this  model  is  the  same  as  that  employed  in  the  core  model, 
the  new  procedure  employs  a  completely  different  computational  flow 
field  solver  from  that  used  in  both  the  core  model  and  the  continuous 
flow  transient  interaction  model.  This  state-of-the-art  computational 
solver  employs  an  implicit,  upwind,  flux-differenced  total  variation 
diminishing  (TVD)  algorithm,  and  is  extremely  efficient  and  capable  of 
providing  highly-detailed,  time-accurate  descriptions  of  the  plasma  and 
field  properties  in  t.e  interaction  region  of  the  solar  wind  with  the 
Earth's  magnetosphere  at  modest  supercomputer  computational  cost 
(approximately  15  CRAY  XMP  CPU  min/case).  The  new  model  allows  the 
complete  time-accurate  determination  of  the  discontinuous  changes  in  the 
bow  shock  and  magnetosheath  flow  field  that  result  from  the  impingement 
of  an  interplanetary  shock  wave  generated  by  solar  flares  or  other 
occurrences  on  a  previously  established  steady  state  terrestrial  bow 
shock/magnetopause  flow  configuration. 

As  an  indication  of  the  capability  of  the  new  discontinuous  flow 
model,  we  provide  results  in  Figure  11  of  an  application  of  the  model. 
Shown  in  that  figure  are  a  sequence  of  plots  displaying  the  predicted 
plasma  pressure  contours  in  the  magnetosheath  region  that  result  from 
the  impingement  of  an  interplanetary  shock  followed  by  Mach  -  12.0  flow 
on  a  previously  established  steady  Mach  -  8.0  flow  pattern  abut  the 
Earth's  magnetopause.  Beginning  with  the  plot  in  the  upper  left  hand 
corner  of  the  figure,  the  internal  shock  patterns  that  are  created  and 
reflected  between  the  bow  shock  and  the  magnetopause  in  the 
magnetosheath  region  as  the  interplanetary  shock  impacts  on  and  then 
moves  downstream  along  the  flanks  of  Earth's  magnetopause  can  be  readily 
discerned.  Notable  in  the  results  is  the  relatively  rapid  establishment 
of  the  new  steady  state  flow  pattern  as  the  Interplanetary  shock  moves 
past  the  magnetopause. 
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gure  11.  Illustration  of  Results  from  the  Discontinuous  Flow 
Transient  Interaction  Submodel:  Plasma  Pressure  Contours  Resulting 
from  the  Global  Interaction  of  an  Interplanetary  Shock  with  an 
Established  Terrestrial  Bow  Shock/Magnetopause  Configuration 
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In  order  to  expedite  development  of  the  discontinuous  flow 
\  transient  interaction  model,  in  these  initial  calculations  with  the 

model  we  have  held  the  location  of  the  magnetopause  fixed.  This 
simplification,  which  can  be  relaxed  simply  by  introducing  time- 
dependent  grid  metrics,  was  adopted  since  the  addition  of  time-dependent 
grid  metrics  in  the  computational  model  would  have  substantially 
complicated  and  extended  the  development  effort  and  time  for  this  model. 
Since  this  unsteady  interaction  model  will  always  require  substantially 
more  computational  resources  than  the  other  axisymmetric  interaction 
models  developed  and  consequently  will  not  be  employed  for  routine 
interaction  calculations,  we  have  chosen  to  defer  the  introduction  of 
time-dependent  metrics  into  the  model  to  a  later  time.  The  major 
Importance  of  the  discontinuous  flow  transient  interaction  model  to  this 
i  effort  are  the  computational  flow  results  that  it  is  able  to  provide. 

These  flow  results  provide  the  basis  for  construction  of  a  rational, 
computationally-simple.  unsteady  flow  interaction  model  that  can  capture 
the  essence  of  the  discontinuous  interaction  process  involving 
impingement  of  an  interplanetary  shock  on  an  established  terrestrial  bow 
shock  and  magnetopause  configuration  but  at  a  computational  efficiency 
at  least  two  orders  of  magnitude  higher.  Based  on  results  from  the 
current  discontinuous  flow  transient  interaction  model,  we  now  have  the 
means  for  estimating  the  extent  of  and  approximating  the  flow  details 
within  the  transient  flow  region  that  separates  the  initial  and  final 
steady  state  bow  shock  and  magnetopause  configurations.  Those  initial 
and  final  steady  state  configurations  and  associated  flow  conditions  can 
be  directly  determined  by  the  core  interaction  model.  Using  this 
rational  modeling  concept,  the  entire  complex  unsteady  shock-on-shock 
interaction  process  can  be  simulated  by  a  relatively  simple  but  accurate 
and  computationally-efficient  interaction  model. 
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In  summary,  the  significant  accomplishments  that  have  been  achieved 
in  this  modeling  effort  include: 


•  Development  of  a  core  interaction  model  in  user-oriented  format 
for  both  quick  response  use  and  also  for  detailed  scientific 
analysis 


•  Development,  verification,  and  application  of  a  series  of 
extended  interaction  models  including: 


*  Fully  unsteady  B  field  predictive  capability  in  core  model 

*  New  3-D  flow  interaction  model  for  predicting  steady  plasma 
properties  about  arbitrary  3-D  magnetopause  geometries 

*  New  aligned  flow  MHD  interaction  model  for  predicting  steady 
full  MHD  solutions  about  arbitrary  axisymmetric  magnetopause 
geometries 

*  New  flow  interaction  model  for  predicting  continuous  plasma 
transient  phenomena 

*  New  flow  interaction  model  for  predicting  discontinuous 
plasma  transient  phenomena 


The  preliminary  development  of  the  terrestrial  interaction  model  is 
now  successfully  complete.  The  interaction  models  developed,  considered 
in  total,  allow  the  quantitative  prediction  of  details  of  the  plasma  and 
field  properties  associated  with  a  variety  of  interaction  phenomena, 
both  steady  and  unsteady,  to  a  degree  of  accuracy  and  routineness  that 
was  previously  unachievable. 

The  most  important  remaining  work  lies  in  the  further  development 
and  enhancement  of  the  modeling  capability  for  accurately  simulating 
strongly  dynamic  processes  associated  with  the  global  solar  wind- 
terrestrial  environment  system,  in  particular  with  regard  to  plasma 
properties.  This  work  should  involve  development  of  both  advanced 
computational  models  capable  of  accurate  dynamic  simulations  as  well  as 
computational ly-simnii^ied,  rational  dynamic  models  which  are  based  upon 
the  advanced  computational  models  but  which  allow  the  accurate 
approximate  simulation  of  unsteady  solar  wind-terrestrial  interaction 
phenomena  at  a  greatly  reduced  computational  cost. 
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Development  of  a  3-D  Gasdynamic  Model  for  Solar  Wind 
Interaction  with  Nonaxi symmetric  Magnetospheres 

S  S  Stahara  and  R  R  Rachiele  (RMA  Aerospace. 
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J  R  Spreiter  (Div.  of  Applied  Mechanics. 
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2891) 

The  gasdynamic  convected  magnetic  field  model  for  predicting 
solar  wind  flow  past  a  planetary  magnetoionopause  obstacle  has 
been  extended  to  three-dimensions  to  apply  to  obstacles  of 
nonaxisymmetric  shape.  Development  of  the  3-D  computational 
model  is  described,  as  well  as  results  of  an  application  of 
the  new  model  to  the  magnetospheres  of  Jupiter  and  Saturn. 
These  magnetospheres  provide  a  strong  test  of  the  new  3-D 
model  since  the  combined  effects  of  rapid  spin,  large  size, 
and  substantial  ring  current  phenomena  are  believed  to  cause 
the  magnetospheres  of  these  planets  to  be  significantly 
broader  near  the  planetary  equatorial  plane  than  near  the 
noon-midnight  polar  meridian  plane.  With  the  aid  of  the  new 
three-dimensional  model,  it  is  now  possible  to  infer  the 
degree  of  flattening  of  the  magnetospheric  cross  sections  from 
a  knowledge  of  the  locations  of  the  low-latitude  magnetopause 
and  bow  shock  crossings.  In  this  paper,  the  computational 
procedures  of  the  new  model  are  described,  and  calculated 
results  are  presented  for  a  number  of  magnetospheres  of 
elliptic  cross  sections  with  values  ranging  from  1  to  2  for 
the  ratio  a/b  of  major  (equatorial)  to  minor  (polar)  axes. 
This  range  is  sufficient  to  include  values  appropriate  to  both 
Jupiter  and  Saturn.  Comparisons  of  the  model  results  with  the 
locations  of  observed  crossings  of  the  magnetopause  and  bow 
shock  directly  provide  an  estimate  of  the  degree  of  equatorial 
broadening  of  the  magnetospheric  cross  section  for  each  of 
these  planets.  For  Jupiter,  the  results  indicate  a  broadening 
to  a/b  ~  1.75,  a  value  consistent  with  previous  estimates 
determined  from  independent  calculations  of  the  three- 
dimensional  magnetosphere  shape  formed  by  adding  to  the 
planetary  dipole  field  the  magnetic  field  of  an  equatorial 
current  sheet  selected  so  as  to  match  the  observed  Jovian 
magnetic  field  near  the  equatorial  plane.  For  Saturn,  a 
similar  comparison  indicates  a  smaller  broadening  to  a/b  - 
1.25.  The  determination  is  less  certain  than  for  Jupiter, 
however,  because  of  the  smaller  amount  and  greater  scatter  of 
the  available  data. 
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The  development  of  a  new  computational  model  for  predicting  in 
detail  the  mass  loading  phenomena  known  to  occur  when  the 
solar  wind  interacts  with  cometary  and  certain  planetary 
ionospheres  is  described.  The  computational  model  employs  two 
separate  but  coupled  flow  solvers  to  determine  the  steady  flow 
field:  an  asymptotic  time-marching  procedure  for  the  subsolar 
to  terminator  region,  and  spatially-marching  procedure  for  the 
post-terminator  region.  The  procedures  are  based  on  current, 
ful ly- impl i ci t  algorithms  and  employ  fitted  discontinuity 
representations  for  both  the  bow  shock  and  inospheric  obstacle 
surfaces.  A  detailed  description  of  the  computational  model 
is  provided,  in  particular,  the  high  degree  of  accuracy  the 
model  is  able  to  provide  in  regions  near  the  obstacle  surface 
is  described  as  well  as  the  the  high  resolution  of  the 
convected  magnetic  field  prediction  near  the  magnetoionopause 
surface.  Mass  loading  due  to  photoionization  has  been 
incorporated  into  the  model,  and  results  of  initial 
applications  of  '  ..e  model  to  cometary  interactions  and  to  the 
mass  loaded  Venus  ionosheath  are  described.  For  the  Venus 
application.  tl  j  model  has  been  used  in  an  inverse  sense  to 
determine  the  mass  loadings  required  to  displace  the  Venusian 
bow  shock  at  the  terminator  to  the  locations  observed  at  solar 
EUV  maximum  and  minimum. 
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Abstract 

The  gasdynamlc  convected  magnetic  field  model 
for  predicting  solar  wind  flow  past  a  planetary 
magnetolonopause  obstacle  has  been  extended  to 
apply  to  obstacles  of  nonaxl symmetric  shape.  The 
need  for  such  an  extension  Is  small  for 
applications  to  the  Earth  and  other  terrestrial 
lanets;  but  Is  significant  for  Jupiter  and 
aturn  where  the  combined  effects  of  rapid  spin, 
large  size  and  substantial  ring  current  phenomena 
are  thought  to  cause  the  magnetosphere  to  be 
significantly  broader  near  the  planetary  equator¬ 
ial  plane  than  In  the  noon-midnight  polar  meridian 
plane.  Confirmation  by  direct  observation  cannot 
yet  be  made,  however,  because  observational  data 
are  available  only  for  conditions  near  the 
equatorial  planes.  Development  of  the  computa¬ 
tional  procedures  of  the  new  three-dimensional 
model  Is  described  and  results  of  an  application 
to  examine  the  suspected  polar  flattening  of  the 
Jupiter  and  Saturn  magnetospheres  are  provided.  A 
quantitative  determination  of  the  degree  of  polar 
flattening  Is  obtained  for  each  of  these  planets 
by  employing  the  new  model  together  with  observa¬ 
tions  of  the  low  latitude  magnetopause  and  bow 
shock  wave  locations.  For  Jupiter,  the  results 
both  confirm  and  better  define  previous 
qualitative  estimates  of  the  flattening  determined 
from  Independent  calculations  of  three-dimensional 
magnetosphere  shapes  formed  by  adding  to  the 
planetary  dipole  field  the  magnetic  field  of  an 
equatorial  current  sheet  selected  so  as  to  match 
the  observed  Jovian  magnetic  field  near  the 
equatorial  plane.  For  Saturn,  the  predicted 
location  of  the  equatorial  bow  shock  agrees  well 
with  the  observations  near  the  nose  of  the  bow 
shock,  but  the  observed  shock  wave  appears  to 
flare  out  substantially  more  than  the  predicted 
shock  wave  at  the  terminator  location.  The  number 
of  these  observations  Is  small,  however,  and  It  Is 
not  possible  at  this  time  to  determine  whether 
this  difference  Is  due  to  some  significant  process 
peculiar  to  Saturn  not  Included  In  the  analysis, 
or  a  meaningless  artifact  resulting  from  temporal 
changes  In  conditions  occurring  at  the  time  the 
observations  were  made. 
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Introduction 

The  gasdynamlc  convected  magnetic  field 
approximation  to  the  magnetohydromagnet 1c  model 
for  supersonic  solar  wind  flow  past  a  planetary 
magnetoionosphere  has  provided  a  useful  basis  for 
understanding  the  solar  wind  planetary  Interaction 
process  from  Its  Inception  more  than  two  decades 
ago*  to  the  present  day  (see  References  2-4  for 
recent  reviews  of  the  theory  and  applications). 
The  original  model  calculations,  although 
accurate,  were  tedious  to  perform  with  the  methods 
available  at  the  time.  As  a  result,  solutions 
were  carried  out  for  only  a  limited  number  of 
cases. 

More  recently,  Sprelter  and  Stahara* •* 
completely  revised  the  solution  procedures  without 
changing  the  underlying  physics  in  order  to  take 
advantage  of  the  enormous  advances  In  computer 
capability  and  computational  methods  that  have 
occurred  since  the  original  model  was  conceived. 
The  result  Is  a  well -documented  model7  capable  of 
routinely  providing  solutions  for  the  detailed 
flow  and  magnetic  field  properties  throughout  the 
magnetolonosheath  region  between  the  bow  wave  and 
a  rather  arbitrarily  specified  shape  for  the 
planetary  magnetopause  or  lonopause.  The  results 
provided  by  the  model  have  now  been  widely  used  in 
the  interpretation  of  spacecraft  data,  and  for 
further  extensions  of  the  theory. 

Throughout  this  theoretical  development  of 
the  model,  the  shape  of  the  magnetospherlc 
obstacle  has  always  been  approximated  as  axl sym¬ 
metric  about  a  line  parallel  to  the  Incident  solar 
wind  direction  and  passing  through  the  planetary 
center.  For  application  to  the  earth,  this  has 
been  justified  by  numerous  studies  of  extensive 
observational  data  In  which  It  has  been  shown  that 
any  departure  from  axlsymmetry  Is  virtually 
Immeasurable.  More  precisely,  Fairfield*  and 
Holzer  and  Slavin’  have  determined  the 
eccentricity  c  of  the  magnetospherlc  cross  section 
In  the  terminator  plane  to  be  less  than  0.2.  Here 
the  eccentricity  Is  defined  as  the  positive  root 
of  c  ■  (1  -  b*/a*  J1'1  where  a  and  b  are  the  major 
and  minor  axes  of  an  ellipse  that  best  fits  the 
observed  magnetospherlc  cross  section  in  the 
terminator  plane.  Observe  that  -c  •  0.2 
corresponds  to  a  value  of  1.02  for  a/b,  t.e.,  the 
approximating  elliptic  cross  section  Is  very 
nearly  circular.  Most  Importantly,  predicted 
results  based  upon  an  axlsymnetrlc  magnetopause 
shape  have  repeatedly  been  found  to  be  In  good 
agreement  with  observations  of  the  location  of  the 
bow  wave  and  of  the  plasma  and  magnetic  field 
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properties  throughout  the  magnetolonosheath  of  the 
earth.  The  same  result  has  also  been  found  at  the 
other  terrestrial  planets." 

Although  only  limited  data,  all  obtained  near 
the  planetary  equatorial  plane,  and  associated 
theoretical  analyses  are  currently  available  for 
Jupiter  and  Saturn,  the  agreement  between  the 
predictions  and  observations  has  been  found  to  be 
clearly  Inferior  to  that  for  the  terrestrial 
planets.  In  a  comprehensive  study  of  the 
positions  of  the  Jupiter  and  Saturn  bow  shocks, 
Slavln  et  « 1 . *  showed  that  the  bow  shock  locations 
predicted  by  the  axlsymmetrlc  model  using  the 
observed  equatorial  shape  of  the  magnetopause  are 
substantially  farther  from  the  magnetosphere  nose 
than  those  actually  observed.  As  shown  In  Figure 
1,  the  resulting  magnetosheath  thickness  In  the 
subsolar  region  Is  45%  too  great  for  Jupiter  and 
20%  too  great  for  Saturn.  Slavln  et  al. 
attributed  these  differences  to  a  substantial 
departure  from  axl symmetry  In  the  magnetospherlc 
shapes  for  these  planets.  In  this  view,  the  polar 
regions  of  the  magnetopause  are  considerably 
flattened  with  respect  to  the  equatorial  regions 
by  effects,  presumably,  of  a  large  equatorial 
current  system"  brought  about  by  the  substantial 
centrifugal  effects  associated  with  the  high  spin 
rates  and  large  size  of  these  planets. 

On  the  basis  of  the  available  Information, 
Slavin  et  al.  suggested  the  polar  flattening 
hierarchy  Illustrated  in  Figure  2  for  the 
magnetopause  cross  sections  in  the  terminator 
plane.  As  shown,  the  cross  section  of  the  earth's 
magnetosphere  is  nearly  circular  with  0.98  <  b/a  ± 
I. 00  corresponding  to  an  eccentricity  e  <  0.2, 
whereas  that  for  the  Jovian  magnetosphere  is 
flattened  into  approximately  an  ellipse  with  b/a  * 
0.6  corresponding  to  c  *  0.8.  The  corresponding 
values  for  Saturn  were  not  determined  to  an  equal 
precision,  but  it  was  concluded  that  this  planet 

?  resents  an  Intermediate  case  with  0.6  <  b/a  < 
.00,  or  c  <  0.8. 

The  polar  flattening  hypothesis  is  consistent 
with  most  of  the  results  shown  In  Figure  I  since 
the  axlsymmetrlc  shape  that  results  from  rotating 
the  observed  equatorial  magentopause  trace  about 
the  symmetry  axis  produces  a  shape  that  has 
significantly  greater  cross-section  area  than  a 
flattened  shape  having  the  same  equatorial 
magnetopause  trace.  As  a  result,  the  bow  shock 
about  an  axlsymmetrlc  magnetopause  will  be 
everywhere  farther  from  the  magnetopause  than  it 
Is  for  a  polar  flattened  magnetopause  having  the 
same  equatorial  shape,  Just  as  observed  for 
Jupiter  and  near  the  magnetopause  nose  for  Saturn. 
The  reason  for  the  pronounced  flaring  of  the 
Saturnian  bow  wave  near  the  terminator  plane  as 
shown  In  Figure  1  Is  not  known,  but  the  amount  of 
data  on  which  this  result  Is  based  Is  so  small 
that  It  may.  In  fact,  be  an  erroneous  steady-state 
interpretation  of  data  obtained  over  times  of 
significant  temporal  variations. 

To  provide  the  appropriate  theoretical 
solutions  for  application  to  such  flattened  magne¬ 
tospheres,  we  have  extended  our  previous  analysis 
so  that  It  can  be  applied  to  the  Interaction  of 
the  solar  wind  flow  with  a  general  three-dimen¬ 
sional  magnetoionospheric  obstacle.  This 


extension  serves  the  dual  purpose  of  providing  (a) 
more  accurate  theoretical  results  required  for 
comparative  studies  with  observations  of  the  solar 
wind  interaction  with  substantially  nonaxisym- 
mclrlc  magnetospheres ,  such  as  thoso  of  Jupiter 
and  Saturn,  and  (b)  an  important  element  of  a 
general  interactive  scheme  In  which  the 
three-dimensional  shape  of  the  magnetoionopause 
and  the  surrounding  flow  and  magnetic  fields  for 
any  planet  are  calculated  alternately  In  the 
process  of  converging  toward  an  exact  solution  of 
the  magnetohydrodynamic  equations. 

In  this  paper,  the  new  computational  model  is 
described  and  results  are  presented  for  the  shape 
of  the  bow  wave  and  the  flow  properties  In 
magnetosheaths  of  polar  flattened  magnetopauses 
representative  of  Jupiter  and  Saturn.  Through  a 
parametric  study  In  which  the  amount  of  polar 
flattening  Is  varied,  a  quantitative  determination 
is  obtained  of  the  degree  of  flattening  at  both 
Jupiter  and  Saturn.  The  new  three-dimensional 
results  are  shown  to  be  In  good  agreement  with  the 
observations,  and  account  for  most  of  the 
differences  between  the  observations  and  the  axi- 
symmetric  model  results. 


Mathematical  Model  and  Solution  Procedure 

The  fundamental  assumption  underlying  the 
present  work,  as  in  the  previous  analyses  of  axi- 
symmetric  flows  reported  in  the  references  cited 
above,  is  that  the  average  bulk  properties  of  the 
solar  wind  flow  around  a  planetary  magnetoiono¬ 
sphere  can  be  described  adequately  by  solutions  of 
the  continuum  equations  of  magnetogasdynamics  of  a 
perfect  gas  having  infinite  electrical 
conductivity  and  zero  viscosity  and  thermal 
conductivity.  The  primary  justification  for  this 
assumption  is  provided  not  by  rigorous  proof  but 
by  the  outstanding  general  agreement  of  the 
calculated  results  with  a  large  amount  of  direct 
observations  in  space. 

Under  this  assumption,  the  equations  used  to 
describe  the  flow  of  the  solar  wind  plasma  past  a 
planetary  magnetoionosphere  are  those  of 
single-fluid  dissipationless  magnetohydrodynamics 
for  the  conservation  of  mass,  momentum,  energy, 
and  magnetic  flux,  augmented  by  the  equation  of 
state  for  a  perfect  gas  and  by  entropy  considera¬ 
tions.  These  equations  have  been  given  many  times 
(see  Reference  2  for  our  most  recent  account),  and 
will  not  be  repeated  here. 

For  typical  solar  wind  conditions,  and 
particularly  at  the  locations  of  the  outer 
planets,  both  the  oncoming  sonic  Mach  number  Ms^  * 
Va./tYP./pos)*/1  and  Alfven  Mach  number  » 
Vt./(B„*/(4*P »))'/*  *re  high  (of  the  order  of  10). 
Under  these  conditions,  an  Important  simplifica¬ 
tion  of  the  magnetohydrodynamic  equations  occurs 
because  the  magnetic  terms  may  be  disregarded  as 
small  in  the  momentum  and  energy  conservation 
equations.  With  these  terms  omitted,  the 
magnetogasdynamic  equations  separate  Into  two 
sets:  the  ordinary  gasdynamics  equations  for  the 
flow  of  a  perfect  gas,  and  the  Faraday  magnetism 
equations.  The  equations  for  the  fluid  motion 
about  an  obstacle  of  given  shape  thereby  reduce  to 
those  of  gasdynamics,  and  can  be  solved  without 


34 


further  consideration  of  the  magnetic  field.  The 
magnetic  field  can  be  determined  subsequently  by 
solving  the  remaining  Faraday  equations  using  the 
values  for  the  velocity  and  density  fields  alroady 
calculated.  The  magnetic  field  determined  In  this 
fashion  Is  usually  Interpreted,  somewhat 
ambiguously,  as  being  'frozen-in'  or  'moving  with 
the  fluid*.  The  resulting  equations  and  the  model 
they  represent  are  commonly  referred  to  as  the 
gasdynamlc  convected  field  approximation.  It 
should  be  noted,  however,  that  the  Interpretation 
of  the  field  as  frozen-in  or  moving  with  the  fluid 
Is  not  dependent  upon  the  separation  of  the  fluid 
motion  determination  from  the  magnetic  field 
determination,  but  Is  actually  an  exact  conse¬ 
quence  of  the  basic  assumption  that  the  flow  can 
be  represented  by  the  magnetohydrodynamic  equa¬ 
tions  for  a  perfect  dissipationless  electrically 
conducting  gas. 

In  carrying  out  the  decoupled  gasdynamlc 
calculation,  the  shape  of  the  magneto ionopause  is 
not  known  a  priori,  but  must  be  found  as  part  of 
the  overall  solution.  In  our  previous  analyses, 
the  magnetoionopause  shape  has  been  determined  by 
either  an  Independent  approximate  theoretical 
analysis,  a  synthesis  of  actual  observational 
results,  or  a  combination  thereof.  This  latter 
procedure  has  also  been  employed  In  the  present 
analysis. 

Computational  methods  based  on  the  full 
magnetohydrodynamic  equation  set  now  exist,  but 
are  still  In  a  preliminary  state  of  development 
and  lack  adequate  resolution  for  many  purposes. 
Hence,  Russell1’  has  recently  stated,  “Currently, 
the  only  way  to  achieve  the  spatial  resolution 
needed  for  a  useful  comparison  with  data  obtained 
on  a  pass  through  a  planetary  magnetosheath  is  to 
employ  the  gasdynamlc  convected  field  model". 

One  of  the  most  Important  and  unique  features 
of  both  the  current  three-dimensional  and  previous 
axlsymmetrlc  models  Is  the  use  of  two  separate, 
but  coupled,  computational  procedures  for 
determining  the  steady  flow  field.  As  Illustrated 
In  Figure  3,  the  model  Is  composed  of  an  obstacle 
nose  region  solver  that  determines  the  flow  field 
from  the  subsolar  region  to  or  somewhat  beyond  the 
temlnator  location  by  using  an  unsteady  procedure 
that  marches  the  solution  forward  In  time  to 
obtain  the  asymptotic  steady  state  solution.  That 
solver  Is  then  linked  together  with  a  spatially- 
marching  tall  region  solver  that  advances  the  flow 
field  solution  downstream  to  any  arbitrary 
distance  In  which  the  flow  remains  supersonic. 

The  precise  axial  location  at  which  the  two 
solutions  are  joined  Is  relatively  arbitrary,  with 
the  single  requirement  being  that  at  all  points  on 
the  joining  plane  the  axial  component  of  the  local 
sonic  Mach  number  be  supersonic.  This  Is  the 
basic  requirement  of  the  spatially-marching 
solver.  In  most  of  our  previous  applications.  It 
was  convenient  to  select  for  the  joining  location 
the  plane  through  the  planetary  center  oriented 
normal  to  the  oncoming  solar  wind  direction,  l.e., 
the  terminator  plane.  In  some  cases  Involving 
magnetoionopause  obstacles  that  possess  consider¬ 
able  flare  at  the  terminator,  it  Is  necessary  to 
locate  the  joining  plane  at  some  fraction  of  an 
obstacle  nose  radii  downstram  of  the  temlnator  In 


order  to  satisfy  the  requirement  of  supersonic 
axial  Mach  number. 

Tho  reasons  for  expanding  the  substantial 
additional  development  effort  to  achieve  such  a 
coupled  two-flow  solver  combination,  rather  than 
employ  a  single  time-marching  solver  for  the 
entire  flow  field  as  is  done  In  all  existing 
magnetohydrodynamic  solutions  procedures,  are 
twofold.  First,  there  Is  the  huge  gain  In 
computational  efficiency  since  the  computationally 
expensive  time-marching  procedure  required  to 
treat  the  subsonlc/transontc  flow  In  the  noise 
region  Is  used  only  where  needed,  while  the 
remainder  of  the  solution  Is  obtained  using  a 
highly  efficient  spatially  marching  procedure 
appropriate  for  supersonic  flow.  Secondly,  the 
computational  capability  is  attained  to  determine 
flow  conditions  as  arbitrarily  far  downstream  as 
needed.  Such  a  capability  cannot  be  provided  by  a 
time-marching  procedure.  Finally,  we  note  that 
such  a  subdivison  of  both  flow  field  and  solution 
procedures  Is  not  uncommon  In  calculating  super¬ 
sonic/hypersonic  flows  and  has  been  used 
extensively  not  only  in  all  our  previous  work,  but 
also  for  example.  In  the  analysis  of  reentry  flows 
including  space  shuttle  applications. 


Nose  Region  Solution  Proc e dure 


In  the  nose  region,  we  have  adapted  the 
procedure  of  Rizk  et  al.  to  solve  the  unsteady, 
three-dimensional  dissipationless  gasdynamlc 
partial  differential  equations  governing  the 
conservation  of  mass,  momentum,  and  energy.  Upon 
introduction  of  the  generalized  independent  vari¬ 
able  transformation  r  »  t,  C  «  CU.y.z),  o  * 
n(x,y,z),  c  “  C ( x ,y , z ) ,  these  equations  can  be 
written  in  strong  conservation  law  form  as: 


(1) 
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Here  (u.v.w)  are  the  Cartesian  velocity 
components,  p  and  p  are  density  and  pressure,  e  is 
the  Internal  energy  and  is  related  to  pressure  by 
p  -  (v-l)le-P(u  ♦  vJ  ♦  w*  )/2] ,  where  y  Is  the 
ratio  of  specific  heats,  (U.V.W)  are  the  contra- 
variant  velocities  without  metric  normalization 
and  are  given  by 


(3) 
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where  (Cx.Cy.C2).  (nx.ny,n»),  Ux.Cy.Cz)  #re  the 
transformation  metrics  and  J  Is  the  Jacobian  of 
the  transformation.  At  each  time  step  the 
trana format Ion  metrics  are  unknown  and  must  bo 
determined  as  part  of  the  solution.  In  our 
Implementation,  they  are  evaluated  numerically 
using  second-order  central  differences.  The  time 
dependent  metric  terms  are  functions  of  the 
Instantaneous  shock  speed,  which  Is  determined 
from  the  shock  fitting  procedure  of  Kutler  et 
al.»- 

The  numerical  algorithm  used  to  solve  the 
above  equations  Is  based  on  the  Beam  and  Warming13 
procedure.  The  method  Is  employed  In  delte  form 
and  Is  fully  Implicit,  noniterative,  second-order 
accurate  In  space,  and  spatially  factored.  When 
applied  to  equation  (1),  the  algorithm  can  be 
written  as: 

[I+AT«tA-CiJ*1VcAtJJ  •  [I+AT6nB-c1J_,VnAnJ]  • 


[I+AT6cC-eiJ‘1VcAcJ]Aq  *  -At[6^E  ♦  6nF  *  i^G] 

-  ceJ*‘[(VcAt)J  ♦  (VnAn)J  +  (VcA;)J]Jq  (4) 

where  (A.B,C)  are  the  Jacobian  matrices  A  =  3E/3q, 
B  ■  3?/3q,  £  ■  3fi/3q,  Aq  *  (q/J)n+1  -  (q/J)n  where 
n  represents  the  time  level,  (6£,6n,6£)  represent 
second-order  central  difference  operators,  (ci.ce) 
are  Implicit  second-order  and  explicit 
fourth-order  smoothing  parameters,  (V, a)  represent 
first-order  backward  and  forward  difference 
operators,  and  t  Is  the  time  step.  A  single  time 
step  integration  of  equation  14)  consists  of  three 
separate  Inversion  steps  associated  with  each  of 
the  three  spatial  directions  (£,n,C).  Each 
Inversion  step  Involves  the  solution  of  a  5x5 
block  tridiagonal  system.  Integration  step  size 
Is  established  by  usinq  the  maximum  eigenvalue  of 
the  Jacobian  matrices  (A,8,C). 

The  analysis  Is  initiated  by  Introducing  a 
computational  mesh  to  discretize  the  flow  field. 
For  the  present  three-dimensional  model  as  well  as 
In  all  our  previous  work,  we  employ  a  fitted 
discontinuity  surface  representation  of  both  the 
bow  shock  and  the  magnetospheric  obstacle.  For 
application  to  three-dimensional  planetary  magne¬ 
tospheres  we  have  found  it  convenient  to  develop  a 
mesh  generator  based  upon  a  generalized  spherical 
(r,0,4)  coordinate  system  with  the  origin  usually 
located  at  the  planetary  center,  but  occasionally 
shifted  downstream  from  that  location  on  a  line 
through  the  planetary  center  and  parallel  to  the 
oncoming  solar  wind. 

The  associated  Cartesian  (x,y,z)  coordinate 
axes  and  a  typical  grid  employed  In  our  analysis 
are  Illustrated  In  Figure  4  where  a  perspective 
cut  out  view  of  one  quadrant  of  the  flow  field  Is 
provided.  In  that  figure  the  origin  coincides 
with  a  location  approximately  half  of  an  obstacle 
radius  downstream  of  the  planetary  center  and  the 
x-axIs,  which  coincides  with  the  longitudinal  axis 
of  the  magnetoionosphere,  points  directly  Into  the 
oncoming  solar  wind.  The  y  axis  lies  In  the 
planetary  equatorial  plane,  and  the  z  axis  points 
In  the  north  polar  direction.  In  Figure  4,  the 


obstacle  to  shock  (r,6)  grids  In  the  azimuthal 
planes  containing  the  x-z  (4  ■  0*)  and  x-y  ($  ■ 
90*)  axes  are  shown  together  with  the  mesh 
distribution  on  the  magnatoionopauta  surface,  and 
the  outflow  boundary  grid  located  In  the  y-z  plane 
at  the  final  downstream  x  location. 

The  quadrant  of  the  (r,6,$)  mesh  shown  in 
Figure  4  contains  (19,30,9)  points,  with  the  9 
azimuthal  planes  spaced  at  equal  10°  Increments. 
In  each  azimuthal  plane,  the  capability  of 
generating  a  body  normal  mesh  has  been  Implemented 
In  the  grid  generator.  A  clustering  capability 
has  also  been  Implemented  to  enable  Independent 
grid  point  clustering  In  each  of  the  three 
coordinate  directions.  For  the  magnetospheric 
obstacle  shown  In  Figure  4,  which  represents  a 
polar  flattened  Jupiter  magnetopause  with  a/b  ■ 
1.75,  the  equal  angular  distribution  of  the 
azimuthal  planes  results  In  a  natural  clustering 
of  the  grid  spacing  on  the  body  In  the  polar 
region  In  contrast  to  that  In  the  equatorial 
region.  This  can  be  retained  or  easily  modified 
If  desired  by  employing  the  clustering  option. 

Boundary  conditions  necessary  for  specifying 
a  properly  posed  mathematical  problem  are  that  the 
flow  satisfy  the  unsteady  Rankine-Hugoniot  shock 
relations  along  the  moving  bow  shock  surface,  be 
entirely  supersonic  at  the  downstream  outflow 
boundary,  be  symmetric  about  the  stagnation 
streamline  at  the  nose,  and  be  tangential  to  the 
obstacle  at  its  surface.  Initial  flow  field 
conditions  are  determined  by  use  of  an 
approximating  surface  for  the  bow  shock  wave,  and 
by  prescribing  a  modified  Newtonian  pressure 
distribution  on  the  surface  of  the  magnetoiono- 
pause.  Since  the  maximum  entropy  streamline  wets 
the  obstacle,  that  fact  plus  the  flow  tangency 
condition  on  the  magnetoionopause  serve  to 
determine  the  remainder  of  the  flow  properties  on 
that  surface.  A  linear  variation  for  the  flow 
properties  between  the  bow  shock  and  the  magneto¬ 
ionopause  then  provides  the  starting  flow  field 
which  is  then  integrated  in  a  time-asymptotic 
fashion  until  the  steady  state  solution  is 
obtained. 

At  the  boundaries,  modification  of  the 
differencing  algorithm  to  account  for  the  appro¬ 
priate  conditions  described  above  is  accomplished 
as  follows.  For  the  fitted  discontinuity  surface 
representing  the  bow  shock,  the  flow  variables 
immediately  downstream  of  the  shock  are  determined 
using  the  unsteady  shock  fitting  procedure  of 
Kutler  et  al.*1'  At  the  outflow  boundary  where  the 
flow  Is  entirely  supersonic,  the  dependent 
variables  are  determined  by  extrapolation  from  the 
adjacent  Interior  points.  Along  the  stagnation 
streamline,  which  forms  a  singularity  line  due  to 
our  use  of  a  spherical-like  coordinate  system13, 
symmetry  conditions  are  enforced  by  employing  a 
second-order  accurate  averaging  procedure.  Also 
Incorporated  In  the  numerical  solver  Is  t ite 
capability  to  enforce  symmetry  boundary  conditions 
across  any  two  arbitrary  azimuthal  planes.  For 
the  present  application,  all  members  of  the  family 
of  three-dimensional  obstacles  being  considered 
possess  quadrilateral  symmetry  about  two 
orthogonal  azimuthal  planes  passing  through  the 
equatorial  and  polar  meridians.  Consequently,  the 
flow  calculations  can  be  restricted  to  a  single 


quadrant  as  shown  In  Figure  4,  thereby  reducing  by 
a  factor  of  four  the  number  of  grid  points  at 
which  the  solution  must  be  determined.  The 
obstacle  surface  flow  tangency  condition  Is 
Incorporated  by  setting  the  contravarlant  velocity 
component ‘normal  to  the  surface  equal  to  zero,  and 
then  det'ermlnlng  the  other  two  contravarlant 
velocity  components  and  the  density  at  the 
obstacle  surface  by  extrapolation  from  Interior 
points.  The  total  energy  is  then  evaluated  by 
using  that  Information  together  with  the  pressure 
computed  from  the  normal  momentum  equation.  The 
Interior  flow'  field  bounded  by  these  various 
boundaries  Is  treated  In  shock-capturing  fashion 
and  therefore  allows  for  the  accurate  determina¬ 
tion  of  secondary  shocks  should  any  occur. 


where 

0  -  D 


(P  -  inb  ♦  n(Rs-Rb>E)0  -  I«br  ♦  n(Rs-Ab>r)2) 


G  -  6 


Tall  Region  Solution  Procedure 


In  the  tall  region,  we  have  employed  the 
shock  capturing  technique  of  Kutler  et  al.w  to 
solve  the  steady,  three-dimensional  dlsslpation- 
less  gasdynamlc  equations.  The  conservation 
equations  for  mass  and  momentum  written  In  weak 
conservation  law  form  In  cylindrical  (x,R,$) 
coordinates  can  be  written  as: 


&x  *  ?R  *  fy,  +  H  «  0 


(5) 


where  the  four-component  vectors  (O.P.fi.fl)  are 
defined  by 


The  numerical  algorithm  chosen  to  solve 
equation  (8)  Is  MacCormack's  explicit  second-order 
accurate  predictor-corrector  method.  As  applied 
to  equation  (9),  the  algorithm  can  be  written  as: 

Ui.j  =  Ui.j  '  +  ‘  *  Sc  (51,j+l 


-  «?J>  * 

=  (4hu?  4  ♦  o!  *  ■ 


V’*  =  (i.){u"  .  ♦  u.  .  -  (F*  .  - 

'2‘  1,j  i  ,j  On  '  1  ,j 


‘  Sc  (®i,j+l}  ’  -  05H*j) 


where 


*  un  (ndC.iAn.Jc),F^j  -  F(u"fJ.nAC.lAn,jAc) 


Fi.J  *  F^i ,j*(n+1^C.10n.jc) .  etc. 


In  which  (u,v,w)  denote  velocity  components  In  the 
(x,R,$)  coordinate  directions.  The  governing  set 
of  equations  Is  completed  with  the  addition  of 
energy  conservation  as  given  by  the  equation  for 
the  total  enthalpy 

Ht  ■  h(p,p)  ♦  qJ/2  «  const  (7) 

Here  q  Is  the  magnitude  of  the  velocity  vector  and 
h(o,p)  represents  the  static  enthalpy  which  for  a 
perfect  gas  Is  given  by  (y/(y-1)J (p/p). 

The  governing  equations  are  transformed  to  a 
computational  space  using  the  Independent  variable 
transformation  £  •  x,  n  ■  (R-Rb)/(R$-Rb).  C  •  ♦, 
where  (Rs,Rb)  represent  the  cylindrical  radial 
coordinate  of  the  bow  shock  and  body  surface, 
respectively,  and  *  Is  the  azimuthal  angle.  This 
yields  the  conservation  equation 

♦  Fn  ♦  Gc  ♦  H  •  °  (8) 


Subsequent  to  each  Integration  step  of 
equation  (10)  with  respect  to  the  hyperbolic 
marching  coordinate  £ ,  the  physical  flow  variables 
(p,0,u,v,w)  must  be  decoded  from  the  components  ui 
of  U.  This  necessitates  the  solution  of  five 
simultaneous  nonlinear  equations  consisting  of 
equation  (7)  together  with  the  four  elements  ut  ■ 
(pu,  p+pu*,  puv,  puw).  For  a  perfect  gas,  the 
procedure  yields  a  quadratic  equation  that  can  be 
solved  for  the  supersonic  velocity  component  u. 
The  remaining  physical  variables  can  then  be 
determined  directly. 

At  the  bow  shock,  the  fitted  shock 
discontinuity  approach  of  Thomas  et  al.,T  Is 
employed,  while  at  the  obstacle  surface  the 
surface  tangency  condition  Isanforcedby  applying 
the  method  of  Abbett.  Integration  step  size  In 
the  marching  direction  Is  redetermined  at  each  C 
location  so  as  to  allow  the  maximum  possible  step 
size  consistent  with  stability. 

The  numerical  calculations  for  which  results 
are  reported  here  were  performed  on  the  NASA  Ames 
Research  Center  CRAY  XMP  computer  facility. 
Approximately  20  to  30  minutes  of  CPU  time  were 
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typically  required  for  the  calculations  for  a 
single  case  of  the  three-dimensional  magnetoiono¬ 
sphere  shapes  considered  In  this  study. 
Essentially  all  of  this  time  Is  spent  with  the 
nose  region  aolvor.  For  tha  family  of 
three-dimensional  shapes  studied,  the  tall  region 
solver  Is  extremely  rapid  and  can  carry  the 
solution  downstream  10  obstacle  nose  radii  In 
typically  less  than  10  CPU  seconds.  As  a  final 
runtime  comparison,  we  note  that  the  typical  CRAY 
XMP  CPU  time  required  for  a  single  case  by  the 
combined  nose  and  tall  axl symmetric  model  to 
determine  an  analogous  flow  field  Is  approximately 
15  seconds. 


Results  and  Applications  to  Jupiter  and  Saturn 

Figure  5  shows  the  calculated  results  for  the 
magnetopause  and  bow  wave  locations  for  the 
conditions  suggested  by  Slavln  et  al.1  as 
representative  for  Jupiter.  These  are  that  the 
flow  has  a  free- stream  sonic  Mach  number  of  10 
and  ratio  of  specific  heats  y  of  2.0,  and  the 
magnetopause  has  an  elliptic  cross  section  with 
constant  ratio  of  major  to  minor  axes  a/b  of  1.75, 
corresponding  to  an  eccentricity  e  of  0.821.  The 
equatorial  magnetopause  profile  Is  taken  to 
approximate  that  deduced  from  observations  at 
Jupiter.  Calculation  of  a  converged  solution  for 
these  conditions  has  been  carried  out  over  the 
region  extending  from  the  nose  of  the  bow  wave  to 
about  1  magnetopause  nose  radius  downstream  of  the 
planetary  center,  or  about  60  Jupiter  radii.  The 
results.  In  the  form  of  traces  of  the  magnetopause 
and  the  bow  shock  wave  In  (a)  each  of  the  ten 
azimuthal  grid  planes  from  the  equatorial  to  the 
polar  planes  of  symmetry  and  (b)  the  terminator 
plane,  display  that  the  flattening  of  the  bow  wave 
Is  much  less  than  that  of  the  magnetopause. 

Details  of  this  three-dimensional  flow  field 
are  presented  In  Figure  6.  In  that  figure, 
contour  plots  of  the  density  P/P«,  pressure  p/p», 
temperature  T/T„,  flow  speed  v/v®,  and  the  sonic 
Mach  number  Ms  In  the  polar  and  equatorial  planes, 
with  an  Isometric  drawing  of  the  magnetopause  and 
bow  shock  wave  In  the  equatorial  and  noon-midnight 
meridian  planes  are  given. 

Figure  7  presents  a  comparison  of  the  Jovian 
magnetopause  and  bow  shock  shapes  observed  near 
the  equatorial  plane1  with  the  bow  shock  shapes 
calculated  for  the  same  equatorial  magnetopause 
shape  for  four  values  for  a/b,  namely  1.00 
(circular  cross  section),  1.50,  1.75,  and  2.00, 
that  bracket  the  value  of  5/3  suggested  by  the 
observations.  The  results  show  that  the  bow  wave 
Is  farthest  out  for  the  axlsymmetrlc  magnetosphere 
and  moves  progressively  nearer  as  the  flattening 
of  the  magnetospherlc  cross  section  Increases. 
This  Is  as  appropriate  since  the  cross-section 
area  of  the  magnetosphere  Is  largest  for  the  axl¬ 
symmetrlc  case,  and  diminishes  progressively  as 
the  magnetopause  Is  flattened  from  a/b  of  1.00  to 
2.00,  thereby  presenting  less  of  an  obstruction  to 
the  flow. 

Comparison  with  the  observations  for  Jupiter 
show  that  they  are  In  excellent  agreement  with  the 
calculated  results  for  a/b  •  1.75,  the  magneto- 
spheric  shape  with  elliptic  cross-section  closest 


to  that  proposed  for  Jupiter  by  Slavln  et  al.1 
This  finding  Is  also  consistent  with  the  result  of 
Engle  and  Beard  In  which  the  three-dimensional 
shape  of  the  Jovian  magnetosphere  was  calculated 
for  a  planotory  dlpOle  fiold  supplemented  by  an 
equatorial  current  sheet  chosen  to  reproduce  the 
observed  Jovian  magnetic  field  near  the  equatorial 
plane.  All  of  this  confirms  the  general 

expectation  that,  provided  the  solution  Is 
determined  for  a  suitably  flattened  magnetopause, 
the  gasdynamlc-convected  field  model  should 
provide  an  even  better  representation  for  the 
magnetosheath  flow  at  Jupiter  than  at  earth 
because  of  the  high  sonic  and  Alfvenlc  Mach 
numbers  that  prevail  there. 

Figure  8  shows  the  corresponding  results  for 
conditions  suggested  by  Slavln  et  al.  as 
representative  of  Saturn.  These  are  that  the  flow 
has  a  free-stream  sonic  Mach  number  Ms-  of  12,  a 
value  of  2.0  for  the  ratio  of  specific  heats  y, 
and  is  about  a  magnetopause  having  the  Indicated 
shape  In  the  equatorial  plane  and  an  elliptic 
cross  section  with  constant  ratio  of  major  to 
minor  axes  a/b  somewhere  between  that  for  the 
earth  and  Jupiter.  Calculated  equatorial  plane 
traces  of  the  bow  shock  wave  are  presented  for 
magnetopause  shapes  with  a/b  •  1.00  (circular 

cross  section),  1.15,  1.25,  1.35,  and  1.45 

together  with  a  best-fit  curve  to  the  observations 
at  Saturn1. 

The  results  for  a/b  *  1.25  match  the  observa¬ 
tions  well  near  the  nose  of  the  bow  shock  wave, 
but  depart  significantly  before  the  planetary 
terminator  plane  is  reached.  The  observed  bow 
wave  flares  out  much  farther  from  the  planet  than 
not  only  the  calculated  results  for  a/b  ■  1.2b, 

but  the  results  for  all  the  calculated  cases, 
including  those  for  an  axlsymmetrlc  magnetopause. 
This  observation  Is  surprising  In  view  of  the 
generally  satisfactory  agreement  between  the 
results  of  the  gasdynamic  calculations  and  the 
observations  for  all  the  other  planets,  and  the 
reason  for  it  remains  unknown.1  It  may  be 
indicative  of  some  significant  physical  process 
not  Included  In  the  analysis,  or  It  might  simply 
be  an  Insignificant  artifact  resulting  from  the 
interpretation  of  a  small  number  of  observations 
and  of  temporal  changes  In  conditions  occurring  at 
the  time  they  were  made. 


Conclusions 

The  development  of  a  computational  model  for 
determining  the  solution  of  the  gasdynamic  portion 
of  the  gasdynamic  convected  magnetic  field  model 
for  solar  wind  flow  past  three-dimensional  magne- 
toionopauses  of  nonaxl symmetric  shape  Is 
described.  Specific  application  of  the  model  is 
made  to  Jupiter  and  Saturn  where,  as  Inferred  from 
observations,  significant  polar  flattening  of  the 
magnetolonpause  cross  sections  exists.  Results 
from  the  new  three-dimensional  model  converge 
continuously  to  those  calculated  previously  for 
axlsymmetrlc  magnetolonspheres  that  were 
determined  using  completely  different  computa¬ 
tional  procedures,  thereby  verifying  the  accuracy 
of  both  methods. 
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Because  the  cost  of  determining  these 
three-dimensional  solutions  Is  a  hundred  times 
greater  than  for  axlsymmetrlc  flows,  and  the 
differences  between  the  results  are  small  for 
•mall  departures  from  axlsyinnotry,  wo  suygost  that 
the  axlsymmetrlc  model  be  used  whenever  the  magne- 
tolonopause  cross  section  Is  nearly  circular,  and 
that  the  new  nonaxlsymmetrlc  model  be  employed 
only  for  those  cases  In  which  substantial 
departures  from  axlsymmetry  exist.  Examples  of 
the  latter  are  Jupiter  and  Saturn,  for  which 
previous  comparisons  between  the  observations  and 
results  calculated  for  axlsymmetrlc  magnetospheres 
were  shown  to  be  In  substantially  poorer  agreement 
with  observations  than  for  the  terrestrial 
planets. 

Similar  comparison  of  results  from  the  new 
model  for  a  parametric  study  of  three-dimensional 
nonaxlsymmetrlc  magnetospheres  having  elliptic 
cross  sections  of  various  eccentricities  represen¬ 
tative  of  Jupiter  and  Saturn  provide  a  new 
estimate  of  the  amount  of  polar  flattening  present 
In  the  magnetospheres  of  these  planets.  It  Is 
found  that  essentially  perfect  agreement  with  the 
observational  data  Is  obtained  for  Jupiter  when  a 
value  of  1.75  Is  used  for  the  ratio  a/b  of  major 
to  minor  axes  of  the  cross  section.  for  Saturn, 
similarly  good  agreement  with  the  observational 
data  Is  found  In  the  nose  region  when  a  value  of 
1.25  Is  used  for  a/b.  Both  of  these  values  are  In 
good  agreement  with  previous  estimates  for  a/b 
based  on  other  considerations. 

Near  the  terminator  region  of  Saturn, 
however,  scanty  observational  evidence  indicates 
that  the  Saturnian  bow  wave  may  flare  out 
considerably  more  than  Indicated  by  any  of  the 
model  predictions,  including  that  for  the  axisym- 
metrlc  case.  It  remains  unknown  whether  this 
difference  is  Indicative  of  some  significant 
process  peculiar  to  Saturn  not  included  in  the 
analysis,  or  a  meaningless  artifact  of  the 
steady-state  Interpretation  of  the  observations 
caused  by  temporal  or  other  changes  In  conditions 
at  the  time  the  observations  were  made  at  that 
planet. 
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a)  Jupiter  bow  shock  comparison 


Fig.  1.  Observations  of  the  Jovian  and  Saturnian 
bow  shock  and  magnetopause  In  the  plane¬ 
tary  equatorial  plane  compared  with  pre¬ 
dictions  (dashed  lines)  of  the  gasdynamlc 
model  with  an  assumed  axlsymmetrlc 
magnetopause 
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Fig.  2.  A  conceptual  representation  of  magneto¬ 
pause  cross  sectional  shape  based  upon 
observations  at  the  earth,  theoretical 
magnetic  field  models  at  Jupiter,  and 
Interpretation  of  gasdynamlc  modeling 
results 


Fig.  3.  Illustration  of  flow  field  subdivision 
and  coupled  flow  solver  combination 
incorporated  in  the  3-0  solar  wind/plane¬ 
tary  Interaction  model 
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Fig.  4.  Perspective  cutout  view  of  one  quadrant 
of  i  typical  computational  grid  employed 
in  the  3-D  model  for  the  current  study 
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Fig.  5.  Predicted  bow  shock  and  corresponding 
magnetopause  locations  In  the  10  equally 
spaced  azimuthal  planes  employed  In  the 
3-D  model  for  a  polar  flattened  Jovian 
magnetopause  having  an  elliptic  cross 
section  with  a/b  ■  1.75  for  M  "10  and  y 
•  2.0 
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b)  Equatorial  meridian  plane 
Fig-  6.  (concluded) 


Fig.  7.  Comparison  of  3-0  model  predictions  with 
observational  data  for  the  Jovian  bow 
shock  In  the  planetary  equatorial  plane 
for  magnetopause  shapes  having  elliptic 
cross  sections  with  different  ratios  of 
major/mlnor  axis  a/b  representing  various 
amounts  of  polar  flattening 


Fig.  8. 


Comparison  of  3-D  model  predictions  with 
observational  data  for  the  Saturnian  bow 
shock  In  the  planetary  equatorial  plane 
for  magnetopause  shapes  having  elliptic 
cross  sections  with  different  ratios  of 
major/minor  axis  a/b  representing  various 
amounts  of  polar  flattening 
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Abstract 

The  gasdynamic  connected  magnetic  field  model  for 
predicting  solar  wind  flows  past  a  planetary  magne- 
toionopause  obstacle  has  been  extended  to  apply  to 
cometary  or  planetary  obstacles  for  which  significant  ion¬ 
isation  of  atmospheric  neutrals  occurs  in  the  surrounding 
solar  wind  flow.  Because  of  the  low  density  of  the  solar 
wind  plasma,  the  presence  of  a  neutral  particle  has  no  ef¬ 
fect  on  the  solar  wind  flow  until  the  moment  it  becomes 
ionized,  whereupon  it  makes  itself  evident  as  a  mass,  mo¬ 
mentum,  and  energy  addition  to  the  flow.  To  account 
for  the  consequences,  our  existing  gasdynamic  convected 
magnetic  field  model  was  extended  to  include  the  effects 
of  such  additions  to  the  flow.  The  need  for  such  an  ex¬ 
tension  is  of  paramount  importance  for  applications  to 
comets  and  can  be  of  major  significance  for  application 
to  planets,  such  as  Venus  and  possibly  Mars  that  have 
sufficiently  weak  intrinsic  magnetic  fields  that  the  effec¬ 
tive  planetary  magnetopause  would  lie  below  the  top  of 
the  sensible  atmosphere. 

A  description  of  the  computational  procedures  devel¬ 
oped  for  the  new  mass  loaded  model  is  provided  and  re¬ 
sults  from  a  series  of  applications  to  both  cometary  and 
planetary  bodies  which  involve  a  large  degree  of  mass 
loading  are  presented.  In  order  to  verify  the  model,  these 
results  have  purposely  been  obtained  with  a  high  degree 
of  convergence  and  spatial  resolution.  The  results  indi¬ 
cate  that  as  the  mass  loading  is  increased,  the  oncoming 
flow  upstream  of  the  bow  wave  becomes  increasingly  com¬ 
pressed,  with  the  bow  shock  wave  weakening  and  moving 
progressively  farther  upstream  and  away  from  the  obsta¬ 
cle.  The  results  are  in  accord  with  general  physical  con¬ 
siderations  regarding  source  flow  addition  to  supersonic 
flows  and  also  with  previous  calculations  of  lesser  reso¬ 
lution  and  accuracy.  Significant  corrections  are  provided 
for  related  results  published  recently  by  others. 

Introduction 

The  gasdynamic  convected  magnetic  field  approxima¬ 
tion  to  the  magnetohydromagnetic  model  for  supersonic 
solar  wind  flow  past  a  planetary  magnetoionosphere  has 
provided  a  useful  basis  for  understanding  the  solar,  wind 
planetary  interaction  process.  From  its  inception  more 
than  two  decades  ago1  to  the  present  day  (see  e.g.  Refs. 
2,3  for  recent  reviews  of  the  theory  and  applications),  the 
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model  has  provided  an  increasingly  useful  tool  for  analys¬ 
ing  a  number  of  phenomena  arising  from  the  interaction 
of  the  solar  wind  with  planetary  obstacles  throughout  the 
solar  system. 

The  original  model  calculations,  although  accurate, 
were  tedious  to  perform  with  the  methods  available  at 
the  time.  As  a  result,  solutions  were  carried  out  for 
only  a  limited  number  of  cases.  More  recently,  Spre- 
iter  and  Stahara4'*  completely  revised  the  solution  proce¬ 
dures  without  changing  the  underlying  physics  in  order  to 
take  advantage  of  the  enormous  advances  in  computer  ca¬ 
pability  and  computational  methods  that  have  occurred 
since  the  original  model  was  conceived.  The  result  is  a 
well-documented  model*  capable  of  routinely  providing 
solutions  for  the  location  of  the  bow  wave  and  the  de¬ 
tailed  flow  and  magnetic  field  properties  throughout  the 
magnetoionosheath  region  between  the  bow  wave  and  a 
rather  arbitrarily  specified  shape  for  the  planetary  mag¬ 
netopause  or  ionopause. 

The  results  provided  by  the  model  have  proven  use¬ 
ful  in  the  interpretation  of  spacecraft  data  for  all  planets 
from  Mercury  to  Jupiter.  The  ability  of  the  numerical 
model  to  match  the  observed  positions  of  the  bow  waves 
and  magnetic  fields  of  the  Earth,  Mars,  and  Mercury  has 
been  remarkably  good.  However,  notable  discrepancies 
were  evident  in  the  initial  comparisons  with  observations 
for  Venus,  Jupiter,  and  Saturn.  The  discrepancies  for  the 
latter  two  giant  planets  have  been  attributed*  to  a  sub¬ 
stantia]  3-D  effect  associated  with  the  flattening  of  the 
planetary  magnetopause  cross  sections  in  the  polar  re¬ 
gions,  presumably  as  a  result  of  the  high  spin  rate  and 
large  size  of  these  planets;  whereas  the  obstacle  cross  sec¬ 
tion  was  always  assumed  to  be  circular  in  the  computa¬ 
tional  model.  To  enable  the  appropriate  comparison  with 
the  gasdynamic  convected  field  model  to  be  made,  Sta- 
hara  et  al.*  extended  the  basic  computational  model  to 
•ppiy  to  an  arbitrary  3-D  magnetospheric  obstacle,  and 
have  demonstrated  that  the  results  from  the  3-D  model 
are  in  good  accord  with  the  observations  for  Jupiter  and 
Saturn. 

For  Venus,  observations  indicate  that  the  planet  has 
no  significant  magnetic  field  and  that  a  virtually  ax- 
isymmetric  obstacle  shape  results  from  solar  wind  flow 
past  the  planetary  ionosphere,  much  as  conceived  in  the 
orginal  theoretical  model4'*'*.  Numerous  comparisons  of 
the  calculated  results  with  those  actually  observed  at 
Venus  indicate,  however,  that  the  observed  bow 
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wave  tend*  on  the  average  to  litre  out  somewhat  more 
on  the  flanks  then  predicted  by  the  model.  Thi*  he* 
been  generally  interpreted  to  be  the  result  of  effect*  of 
solar  wind  pickup  of  newly  ionised  atmospheric  particle* 
not  included  in  the  model  (aee  Breus”  for  *  recent  re¬ 
view).  Breus”  end  Belotserkovskii  et  al.14  have  recently 
published  result*  of  a  calculation  to  illustrate  the  magni¬ 
tude  of  these  effects,  but  their  results  differ  substantially 
from  those  provided  by  our  model  for  the  same  condi¬ 
tions.  While  their  results  match  the  observations  reason¬ 
ably  well,  we  believe  their  computed  results  are  lacking 
in  accuracy  and  greatly  overestimate  the  effects  of  ion 
pickup. 

Similar,  but  much  more  extensive,  ion  pickup  processes 
become  of  dominant  importance  for  comets,  and,  as  em¬ 
phasised  by  Breus  et  al.”,  much  is  to  be  learned  by  stud¬ 
ies  of  their  similarities  and  differences.  Study  of  the  na¬ 
ture  of  comets  extends  from  antiquity  to  the  present  day 
(see  Ref.  16  for  a  comprehensive  account),  and  has  pro¬ 
vided  many  insights  into  the  properties  of  both  comets 
and  the  interplanetary  environment.  Indeed,  the  first 
suggestion  for  the  continual  existence  of  a  high  speed 
outflow  of  gas  from  the  sun,  the  solar  wind,  was  made  by 
Biermann”'”  on  the  basis  of  his  studies  of  comet  tails. 
Among  the  first  systematic  studies  of  solar  wind  flow  past 
a  comet  in  the  modern  context  is  that  of  Ref.  19  in  which 
the  solar  wind  flow  was  represented  by  the  conservation 
equations  of  gasdynamics  with  additional  source  terms 
included  to  represent  effects  of  mass,  momentum,  and 
energy  addtion  in  the  flow.  These  are  essentially  the 
same  equations  employed  in  the  model  developed  here, 
although  the  expressions  for  the  source  terms  are  quite 
different.  However,  the  present  numerical  solutions  are 
of  significantly  higher  accuracy,  due  largely  to  advances 
in  both  numerical  methods  and  computational  capability 
made  during  the  intervening  twenty  years. 

Schmidt  and  Wegmann*9  subsequently  rephrased  the 
analysis  in  terms  of  the  equations  of  magnetohydro¬ 
dynamics,  must  as  has  been  done  for  the  planetary 
applications1  and  presented  an  inte|e*tiag  but  limited 
set  of  low  resolution  result*.  Because  .he  computational 
cost  of  determining  solutions  of  the  magnetohydrody¬ 
namic  equations  is  considerably  greater  than  it  is  for 
the  gasdynamic  equations,  and  because  the  latter  have 
a  long  history  of  providing  results  of  useful  accuracy  for 
related  applications,  we  have  based  our  model  on  the  gas- 
dynamic  equations  but  have  employed  superior  methods 
for  their  solution.  In  this  way,  our  model  provides  high 
resolution  results  with  almost  minimal  computing  costs 
using  a  modern  supercomputer.  Although  it  is  highly 
unlikely  that  the  flow  about  an  actual  comet  is  either 
steady  or  possesses  any  simple  symmetry,  the  present  ini¬ 
tial  cometary  application  is  based  on  the  assumption  that 
the  flow  is  steady  and  axisymmetric  under  the  hopeful  as¬ 
sumption  that  such  results  will  be  generally  informative 
about  the  basic  nature  of  cometary  flow  fields.  The  bene¬ 
fits  are  at  least  a  hundredfold  reduction  in  the  amount  of 
computing  required  for  a  solution,  and  the  avoidance  of 
a  requirement  to  specify  with  realism  the  unsteady  and 
nonspherical  symmetry  properties  of  the  outflow  of  neu¬ 


tral  gas  from  the  cometary  nucleus,  which  are  in  general 
unknown. 

In  this  paper,  the  new  computational  model  is  de¬ 
scribed  and  result*  are  presented  for  the  shape  of  the 
bow  wave  and  the  associated  properties  for  solar  wind 
flow  past  both  planetary  and  cometary  bodies  exuding 
neutral  particles  which  subsequently  become  ionised.  A 
quantitative  determination  of  the  effects  of  such  ion  addi¬ 
tion  to  the  solar  wind  is  provided  by  presenting  the  results 
of  a  parametric  study  in  which  the  amount  of  ion  addi¬ 
tion  is  varied.  The  new  results  provide  considerable  im¬ 
provement  over  previous  qualitative  and  semi-quantitive 
determinations  of  the  properties  of  solar  wind  flow  past 
a  comet.  For  Venus  they  also  show  that  a  much  greater 
addition  of  new  ions  than  assumed  by  Breus  et  al.1*  and 
Belotserkovskii  et  al.*4  is  required  to  achieve  good  quan¬ 
titative  agreement  with  observations. 

Mathematical  Model  and  Solution  Procedure 

The  fundamental  assumption  underlying  the  present 
work,  as  in  the  previous  analyses  in  the  references  dted 
above,  is  that  the  average  bulk  properties  of  the  solar 
wind  flow  around  a  planetary  magnetoionosphere  can  be 
described  adequately  by  solutions  of  the  continuum  equa¬ 
tions  of  magnetogaadynamics  for  a  perfect  dissipationless 
plasma  having  infinite  electrical  conductivity  and  sero 
viscosity  and  thermal  conductivity.  The  primary  jus¬ 
tification  for  this  assumption  is  provided  not  by  rigor¬ 
ous  proof  but  by  the  outstanding  general  agreement  of 
the  results  predicted  on  this  basis  compared  with  a  large 
number  of  direct  observations  in  space.  These  governing 
equations,  which  have  been  given  previously  (see  Ref.  12 
for  our  most  recent  account),  will  not  be  repeated  here. 
For  application  to  flows  with  significant  pickup  of  newly 
created  ions,  the  conservation  equations  used  to  describe 
the  solar  wind  flow  must  be  augmented  to  include  term* 
to  represent  the  effective  mass,  momentum,  and  energy 
additions  to  the  flow. 

For  typical  solar  wind  conditions,  both  the  oncoming 
sonic  Mach  number  Msoc  =  »oo /y/iv/p  and  Alfven  Mach 
number  Maoo  —  t>oe/\/®Io7(**£ooj  are  substantially 
greater  than  unity.  Under  these  conditions,  an  impor¬ 
tant  simplification  of  the  magnetohydrodynamic  equa¬ 
tions  may  be  introduced  because  of  the  relative  small¬ 
ness  of  the  magnetic  terms  in  the  momentum  and  energy 
conservation  equations.  With  these  terms  omitted,  the 
magnetogasdynamic  equations  separate  into  two  sets:  the 
gasdynamics  equations  for  the  flow  of  a  perfect  gas  with 
mass,  momentum,  and  energy  addition,  and  the  Faraday 
magnetism  equations.  The  equations  for  the  fluid  motion 
about  an  obstacle  of  given  shape  thereby  reduce  to  those 
of  gasdynamics,  and  can  be  solved  without  further  con¬ 
sideration  of  the  magnetic  field.  The  magnetic  field  can 
then  be  determined  subsequently  by  solving  the  remain¬ 
ing  Faraday  equations  using  the  values  for  the  velocity 
and  density  fields  already  calculated.  The  magnetic  field 
determined  in  this  fashion  is  usually  interpreted,  some¬ 
what  ambiguously,  as  being  ’frosen-in’  or  ’moving  with 
the  fluid’.  The  resulting  equations  and  the  model  they 
represent  are  referred  to  a*  the  gasdynamic- corrected 
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field  approximation.  It  should  be  noted,  however,  that 
the  interpretation  of  the  field  a<  frosen-in  or  moving  with 
the  fluid  it  not  dependent  upon  the  separation  of  the  fluid 
motion  determination  from  the  magnetic  field  determina¬ 
tion,  but  it  an  exact  consequence  of  the  basic  astumption 
that  the  flow  can  be  represented  by  the  magnetohydrody- 
namic  equations  for  a  perfect  dissipationlett  electrically 
conducting  gat. 


In  carrying  out  the  decoupled  gasdynamic  calculation, 
the  shape  of  the  magnetoionopause  is  not  know  apriori, 
but  mutt  be  found  as  part  of  the  overall  solution.  In 
our  previous  analyses,  the  magnetoionopause  shape  has 
been  determined  by  either  an  independent  approximate 
theoretical  analysis,  a  synthesis  of  actual  observational 
results,  or  a  combination  thereof.  This  latter  procedure 
has  also  been  employed  in  the  present  analysis. 

Computational  models  based  on  the  full  magnetohy¬ 
drodynamic  equation  set  now  exist,  but  are  still  in  a  pre¬ 
liminary  state  of  development  and  lack  adequate  reso¬ 
lution  for  many  purposes.  Hence,  Russell1  has  recently 
stated,  “Currently,  the  only  way  to  achieve  the  spatial 
resolution  needed  for  a  useful  comparison  with  data  ob¬ 
tained  on  a  pass  through  a  planetary  magnetosheath  is 
to  employ  the  gasdynamic-convected  field  model". 

To  be  more  precise,  the  conservation  equations  for  un¬ 
steady  axisymmetric  gasdynamic  flow  with  mass,  momen¬ 
tum,  and  energy  addition  may  be  written  as  follows: 


8Q  dE  dF 
dt  +  dz  +  dR 


5  +  ff 


(1) 


in  which 
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{  pu  \ 

+  P 

puv 

V(e  +  pW 


f  5mass 

§  _  ^x-momentum 
^R-momentum 
V  ^energy 

where  Q  is  the  conserved  quantity,  E,  F,  and  B  are  the 
corresponding  fluxes,  and  5  represents  the  source  terms 
as  indicated.  The  other  quantities  have  their  usual  mean¬ 
ings,  i.e.,  p  is  the  density,  a  and  v  are  the  velocity  com¬ 
ponents  parallel  to  the  cylindrical  t  and  R  axes,  p  is  the 
pressure,  and  e  is  the  total  energy,  which  is  related  to  the 
other  quantities  by  p  =  (7  —  l)(e  —  p(u*  +  v*)/2)  where 
7  is  the  ratio  of  specific  heats. 


The  remaining  Faraday  equations  for  the  magnetic  field 


S  are: 

M  =  eurl(V  x  B),  divff  =  0  (3) 

where  V  and  B  are  the  velocity  and  magnetic  field  vectors 
respectively. 

In  our  realisation  of  the  solution,  the  discretized  gasdy¬ 
namic  equations  are  solved  on  a  segmented  inner/outer 
curvilinear  grid  system  in  order  to  accurately  account  for 
mass  loading  effects  upstream  of  the  bow  shock.  As  illus¬ 
trated  in  Figure  1,  the  inner  grid  segment  extends  from 
the  obstacle  surface  to  the  fitted  bow  shock,  while  the 
outer  segment  starts  at  the  bow  shock  and  extends  to 
a  location  sufficiently  far  upstream  where  the  mass,  mo¬ 
mentum,  and  energy  additions  due  to  the  source  terms 
are  no  longer  significant.  As  is  evident  from  the  grid 
displayed  in  Figure  1,  we  employ  a  fitted  discontinuity 
procedure  to  represent  both  the  bow  shock  and  the  ob¬ 
stacle  surface.  Although  this  is  far  more  difficult  to  im¬ 
plement  than  simply  allowing  the  solution  algorithm  to 
capture  these  discontinuities,  employing  a  fitted  discon¬ 
tinuity  representation  allows  the  exact  satisfaction  of  the 
shock  jump  and  obstacle  surface  boundary  conditions  at 
surfaces  that  are  precisely  defined,  in  constrast  to  the 
smeared  representation  that  inevitably  results  from  even 
the  best  capturing  methods. 

One  of  the  most  important  features  of  both  the  present 
and  all  of  our  previous  models  is  the  use  of  two  separate, 
but  coupled,  computational  procedures  for  determining 
the  steady  flow  field.  As  illustrated  in  Figures  1  and 
2,  the  model  is  composed  of  a  nose  region  solver  that 
determines  the  flow  field  from  the  subsolar  region  to  or 
somewhat  beyond  the  terminator  location  by  using  an 
unsteady  procedure  that  marches  the  solution  forward  in 
time  to  obtain  the  asymptotic  steady  state  solution.  That 
solver  is  then  joined  with  a  spatial-marching  tail  region 
solver  that  advances  the  flow  field  solution  downstream 
to  any  arbitrary  distance  in  which  the  flow  remains  su¬ 
personic. 

The  precise  axial  location  at  which  the  two  solutions 
are  joined  is  relatively  arbitrary,  the  single  constraint  be¬ 
ing  that  at  all  points  on  the  joining  plane  the  axial  compo¬ 
nent  of  the  local  sonic  Mach  number  be  supersonic.  This 
is  the  basic  requirement  of  the  spatial-marching  solver  in 
order  that  it  be  able  to  march  the  solution  downstream. 
In  most  of  our  previous  applications,  it  has  been  conve¬ 
nient  to  select  for  the  joining  location  the  plane  through 
the  planetary  center  oriented  normal  to  the  oncoming  so¬ 
lar  wind  direction,  i.e.,  the  terminator  plane.  In  some 
cases  involving  obstacles  that  possess  considerable  flare 
at  the  terminator,  it  is  necessary  to  relocate  the  joining 
plane  at  some  fraction  of  an  obstacle  nose  radii  down¬ 
stream  of  the  terminator  in  order  to  satisfy  the  require¬ 
ment  of  supersonic  axial  Mach  number. 

The  primary  reasons  for  expending  the  substantial 
additional  development  effort  to  realise  such  a  coupled 
two  flow  solver  combination,  rather  than  employ  a  sin- 
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gle  time- marching  solver  for  the  entire  flow  field  as  ia 
done  in  all  existing  magnetohydrodynamic  eolutione  pro¬ 
cedure*,  are  twofold.  Fir*t,  there  i*  the  substantial  gain 
in  computational  efficiency  since  the  computationally  ex¬ 
pensive  time-marching  procedure  required  to  treat  the 
subsonic/ transonic  flow  in  the  nose  region  is  used  only 
where  needed,  while  the  remainder  of  the  solution  is  ob¬ 
tained  using  a  highly  efficient  spatially  marching  proce¬ 
dure  appropriate  for  supersonic  flow.  Secondly,  the  com¬ 
putational  capability  is  attained  to  determine  flow  condi¬ 
tions  as  far  downstream  as  desired.  Such  a  capability  can¬ 
not  be  provided  by  a  time-marching  procedure.  Finally, 
we  note  that  such  a  subdivision  of  both  flow  field  and 
solution  procedures  is  not  uncommon  in  calculating  su¬ 
personic/hypersonic  flows  and  has  been  used  extensively 
not  only  in  all  our  previous  work,  but  also  for  example, 
in  the  analysis  of  reentry  flows  including  space  shuttle 
applications. 


Nose  Region  Solution  Procedure 


In  the  nose  region,  we  have  employed  the  procedure  of 
Beam  and  Warming  30  to  solve  equations  (1)  and  (2)  with 
specified  mass,  momentum,  and  energy  addition.  First, 
the  governing  equations  are  transformed  from  the  (*,  A,  i) 
coordinates  to  the  generalised  set  of  coordinates,  ((,  q,  r ). 
This  new  set  of  equations  can  be  written  as: 


8Q  dE  8F 

dr  +  di  +  df, 


=  5+  H 


(4) 


The  vectors  Q ,  £,  and  F  represent  the  transformed  de¬ 
pendent  variables,  and  their  flux  vectors  as: 


=  £  =  i((.E  +  UF),  F=j{r,tE  +  T,RF) 

(5) 

The  transformed  source  terms  are  simply  written  as: 

B  =  jB,  S=jS  (6) 

In  the  above  equations,  ((«,(/», q«,qi i)  are  the  transfor¬ 
mation  metrics  and  J  is  the  Jacobian  of  the  transforma¬ 
tion.  Even  though  the  grid  points  are  free  to  move  in 
time,  the  unsteady  metric  terms  corresponding  to  this 
grid  movement  have  been  neglected  for  the  interior  nodal 
points  since  only  a  steady  state  solution  is  desired.  This  is 
easily  justified  since  the  grid  speeds  vanish  as  the  steady 
state  solution  evolves. 


The  numerical  algorithm  used  to  solve  the  above  equa¬ 
tions  is  based  on  the  Beam  and  Warming11  procedure 
with  the  boundary  conditions  of  Chakravarthy11.  The 
method  is  a  fully  implicit,  noniterative,  scheme  that  is 
second-order  accurate  scheme  in  space.  It  can  be  written 
ia  delta  form(AQ  *  <5n+1  -  <J")  as: 

[/  +  At*(>i-+AtflA-]”A6  > 

1  1  n  (7) 

-At  (S(E  +  6nF  -  B  -  $)"  +  D*(Q) 

Here,  A,  and  B  are  the  standard  fluid  dynamic  Jaco- 
bians,  A  =  8E/d<$  and  6  =  dP/8Q  and  I  is  the  identity 


matrix.  The  superscript  n  implies  evaluation  at  the  nth 
time  level.  Second-order,  central-differences  have  been 
used  for  all  discretized  derivatives  and  a  fourth  order  ar¬ 
tificial  dissipation  term,  Z>,  has  been  added  for  stability. 
The  integration  step  site,  At,  is  generaly  specified  in  the 
begining  of  the  calculation  and  held  constant  for  the  en¬ 
tire  computation. 

As  noted  above,  the  computations  are  carried  out  on  a 
mesh  that  conforms  to  the  fixed  shape  of  the  ob¬ 

stacle  and  adapts  itself  at  each  time  step  of  the  solution  to 
conform  to  the  time- varying  shape  of  the  bow  wave.  The 
associated  cylindrical  (*,  A)  coordinate  axes  and  a  typi¬ 
cal  grid  employed  in  our  analysis  are  illustrated  in  Figure 
2.  The  origin  coincides  with  the  center  of  the  planet  or 
the  cometary  nucleus,  the  x-axis,  which  coincides  with 
the  longitudinal  axis  of  symmetry  of  the  magnetoiono¬ 
sphere,  points  directly  into  the  oncoming  solar  wind,  and 
the  R  axis  extends  from  the  origin  perpendicular  to  the 
x  axis.  The  obstacle  to  shock  grid  lines  and  their  contin¬ 
uations  extend  in  the  radial  r  direction  from  the  origin. 
The  azimuthal  grid  lines  extend  from  the  *  axis  to  the 
A  axis  at  the  terminator  plane  in  such  a  way  that  one 
of  the  grid  lines  conforms  to  the  obstacle  surface  and  an¬ 
other  conforms  to  the  bow  wave.  Various  numbers  of  grid 
lines  may  be  used  for  the  calculations,  but  the  situation 
portrayed  in  Figure  2  in  which  a  finer  mesh  is  employed 
between  the  obstacle  and  the  bow  wave  than  in  the  region 
beyond  the  bow  wave  may  be  taken  as  representative.  In 
this  example,  there  are  35  azimuthal  grid  lines  and  23 
radial  grid  lines  for  the  nose  region.  Values  for  all  flow 
quantities  are  calculated  at  every  point  of  intersection  of 
these  grid  lines. 

The  boundary  conditions  applied  to  the  set  of  govern¬ 
ing  equations  are  outlined  in  Figure  3.  At  the  outer¬ 
most  boundary,  the  flow  is  uniform  and  parallel  to  the  x- 
direction  with  free-stream  values  for  the  density,  velocity, 
and  Mach  number  corresponding  to  those  of  the  undis¬ 
turbed  solar  wind.  At  the  innermost  boundary,  the  flow 
is  tangential  to  the  obstacle  surface.  In  addition,  the  so¬ 
lution  must  satisfy  the  unsteady  Rankine-Hugoniot  shock 
relations  at  the  moving  bow  shock  surface,  be  entirely  su¬ 
personic  at  the  downstream  outflow  boundary  of  the  nose 
region,  and  be  symmetric  about  the  stagnation  stream¬ 
line  upstream  from  the  nose.  Initial  flow  field  conditions 
are  determined  by  use  of  an  approximating  surface  for 
the  bow  shock  wave,  and  by  prescribing  a  modified  New¬ 
tonian  pressure  distribution  on  the  surface  of  the  mag- 
netoionopause.  Since  the  maximum  entropy  streamline 
wets  the  obstacle,  that  fact  plus  the  flow  tangency  con¬ 
dition  on  the  magnetoionopause  serve  to  determine  the 
remainder  of  the  flow  properties  on  that  surface.  A  lin¬ 
ear  variation  for  the  flow  properties  betwen  the  bow  shock 
and  the  magnetoionopause  then  provides  the  starting  flow 
field  which  is  then  integrated  in  a  time-asymptotic  fash¬ 
ion  until  the  steady  state  solution  is  obtained. 

At  the  boundaries,  modification  of  the  differencing  al¬ 
gorithm  to  account  for  the  particular  conditions  described 
above  is  accomplished  as  follows.  First,  an  eigenvalue 
analysis  reveals  the  number  of  boundary  conditions  need- 
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inf  to  be  specified  it  each  boundary.  The  number  of 
boundary  conditions  at  each  boundary  must  be  equal  to 
the  number  of  characteristic  waves  entering  the  compu¬ 
tational  domain  through  that  boundary.  For  the  super¬ 
sonic  outflow  condition,  since  all  eigenvalues  are  asso¬ 
ciated  with  waves  exiting  the  computational  domain  at 
that  boundary,  no  additional  boundary  conditions  need 
be  specified  and  the  governing  equations  are  discretised 
using  a  one  sided  difference  into  the  computational  do¬ 
main  from  this  boundary.  These  equations  are  then  made 
implicit  by  a  simple  linearisation  process. 


Boundary  conditions  other  than  those  associated  with 
supersonic  flow  are  not  as  straightforward.  A  character¬ 
istic  equation  is  specified  for  each  outward  traveling  wave 
and  boundary  conditions  are  specified  until  the  number 
of  equations  matches  the  number  of  dependent  variables. 
These  characteristic  equations  are  a  direct  result  of  the 
eigenvalue  analysis  and  a  further  discussion  of  their  origin 
can  be  found  in  reference  22.  For  the  surface  tangeney 
boundary  condition,  three  characteristic  equations  and 
one  boundary  condition  are  used  since  only  one  eigen¬ 
value  is  associated  with  an  inward  traveling  wave.  The 
boundary  condition  used  in  this  case  is  the  specification 
a  sero  surface-normal  velocity. 


For  the  bow  shock  boundary  condition,  two  additional 
quantities  are  added  to  the  set  of  dependent  variables. 
They  are  zT  and  yT,  the  shock  propagation  velocities 
speeds  in  the  x  and  y  direction  respectively.  An  eigen¬ 
value  analysis  reveals  that  only  one  wave  propegates  from 
the  interior  flow  field  to  the  bow  shock  resulting  in  only 
one  valid  characteristic  equation.  Five  boundary  con¬ 
ditions  must  then  be  specified.  They  consist  of  four 
Rankine-Hugoniot  relations  across  the  bow  shock  and  one 
grid  point  restriction  equation.  The  Rankine-Hugoniot 
relations  are  well  known  and  the  grid  point  restriction 
equation  restricts  the  shock  grid  points  to  lie  only  on 
rays  eminating  from  the  obstacle  center. 


Tail  Region  Solution  Procedure 


If  the  velocity  of  the  flow  in  the  streamwise  coordinate 
direction  ia  greater  than  the  speed  of  sound,  a  method  ean 
be  employed  that  solves  for  the  flowfield  directly  without 
a  costly  time  iteration.  The  method  is  baaed  on  that  of 
Shiff  and  Steger11  and  marches  along  the  streamwise  co¬ 
ordinate  rather  than  time  to  obtain  a  solution.  Since  the 
flow  in  the  tail  region  fulfills  this  criteria,  such  a  method 
is  used.  The  governing  equations  used  here  are  identical 
to  those  of  equations  4,  with  the  exclusion  of  any  time 
derivatives. 


BE  BF 
6(  +  Bf, 


=  $  +  H 


(8) 


If  a  first  order  appoximation  to  the  streamwise  derivative 
is  made  and  a  spatial  linearisation  of  all  quantities  eval¬ 
uated  at  the  new  spatial  location  is  done,  the  numerical 
scheme  can  be  written  in  delta  form  (AQ  =  $,+  ]  -  <$,) 


[M  +  -ti  +  *,/-  _  _  $) 

(9) 

The  subscript  i  denotes  the  axial  plane,  and  the  primed 
quantities  signify  that  the  metric  quantities  have 
frosen  in  the  linearisation.  Again  the  {-operator  implies 
central  differencing  and  a  fourth  order  artificial  dissipa¬ 
tion  has  been  used.  This  tail  region  solver  has  been  built 
in  such  a  way  to  be  as  consistent  as  possible  with  the 
nose  region  solver.  Even  boundary  conditions  are  sim¬ 
ilar  with  the  exception  of  spatial  rather  that  temporal 
linearisation. 


Applications  to  Comets 

To  complete  the  specification  of  the  mathematical 
problem  to  represent  solar  wind  flow  past  a  comet,  it  is 
necessary  to  supply  expressions  for  the  mass,  momentum, 
and  energy  source  terms.  For  the  applications  reported 
here,  we  have  employed  the  model  suggested  by  Schmidt 
and  Wegmann30  and  others  in  which 

{  Smass  ' 

g  _  ^x-momentum 
^R- momentum 
V  ^energy  / 

where 

Sn  ~  4*Wr> e*P(~w) 

and  where  G  is  the  total  number  of  neutral  particles  emit¬ 
ted  from  the  comes  per  second,  a  is  the  photoionization 
rate,  W  is  the  radial  expansion  velocity  of  the  neutral 
particles,  r  is  the  three-dimensional  radial  distance  from 
the  origin,  and  mt  is  the  neutral  particle  mass.  The  phys¬ 
ical  considerations  underlying  the  above  formulation  for 
the  source  terms  are  that  the  ions  are  formed  at  a  con¬ 
stant  rate  9  by  photoionisation  of  the  originally  neutral 
cometary  particles  which  have  been  evaporated  from  the 
nucleus.  These  neutral  molecules  are  all  considered  to 
have  the  same  mass  mc  and  to  move  radially  away  from 
the  comet  with  a  constant  speed  W.  The  electrons  are 
treated  as  a  cool  component  whose  effects  on  the  proper¬ 
ties  of  the  flow  are  negligible. 

With  regard  to  the  computational  considerations  that 
these  source  terms  have  on  the  numerical  method,  we 
note  that  with  G,  s,  Me,  and  W  considered  to  be  con¬ 
stants,  a  typical  conservation  equation  for,  say,  the  con¬ 
tinuity  of  mass  then  becomes  of  the  form 

L.H.S.  continuity  equation  =  Sm.f(r)  (12) 

where  Sm  =  Constant  and  /(r)  =  0(1).  We  note  that 
representative  values  for  the  normalized  constant  $m  for 
cometary  interactions  range  from  as  small  as  0.1  to  as 
large  as  1000  and  even  greater.  The  larger  values  result 
in  an  extremely  stiff  set  of  partial  differetial  equations  and 


(10) 
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poses  ft  ntck  and  ckalleafini  tot  of  any  compalfttio&ftl 
procedure. 

We  have  applied  the  above  computational  model  to  ft 
pftrftmetric  series  of  mass  loadings  for  the  above  cometary 
source  model  with  good  success.  Figure  4  displays  results 
for  the  density  contours  for  a  complete  nose  and  tail  Sow 
field  for  one  such  calculation  employing  cometary  source 
parameters  representative  of  those  suggested  by  Schmidt 
and  Wegmann".  .  That  result  corresponds  to  a  cometary 
mass  loading  with  G  =  6  x  10*'  particles  per  second,  a 
photoionixation  rate  a  of  10~*  per  second,  a  radial  ex¬ 
pansion  velocity  of  1  km/sec,  a  cometary  particle  mass 
of  23.3mp  where  m,  represents  the  proton  mass,  a  so¬ 
lar  wind  density  pm  =  5 m,/em'  and  velocity  =  350 
km/sec,  sonic  Mach  number  Af;K  =  6.0,  and  ratio  of 
specific  heats  7  =  5/3.  The  assumed  cometary  tangen¬ 
tial  discontinuity  surface  was  approximated  as  a  hemi¬ 
sphere/cylinder.  For  perspective  on  this  choice  for  G, 
recent  observations  provide  a  value  of  about  4  x  10"  par¬ 
ticles/sec  for  Giacobinin-Zinner,  a  middle-size  comet" 
and  about  1  x  10"  particles/sec  for  Halley,  a  large  comet. 
The  results  for  density  contour  lines  upstream  of  the  bow, 
indicating  that  considerable  compression  of  the  flow  oc¬ 
curs  upstream  of  the  bow  wave.  Overall,  however,  the 
results  are  quite  similar  to  previously  published  results 
for  planetary  flows  with  no  ion  pickup. 

Figure  5  shows  corresponding  results  for  the  nose  re¬ 
gion  for  four  different  values  for  G  ranging  from  1.4  x  10" 
to  a  substantially  larger  value  of  1.12  x  10".  Also  shown 
on  those  figures  are  the  corresponding  values  of  the  nor¬ 
malized  mass  source  constant  $m  which  range  from  1.2  to 
9.2.  The  plot  for  the  smallest  G  indicates  only  one  density 
contour  line  upstream  of  the  bow  wave,  and  a  compar¬ 
ison  reveals  that  the  results  in  their  entirety  differ  only 
slightly  from  those  without  ion  pickup  about  a  similar 
obstacle.  The  plots  for  larger  G  show  increased  compres¬ 
sion  upstream  of  the  bow  wave.  They  also  show  that  the 
entire  bow  wave  moves  away  from  the  cometary  obstacle, 
greatly  enlarging  the  region  between  the  obstacle  and  the 
bow  wave.  More  specifically,  the  nose  of  the  bow  wave  is 
shown  to  move  from  about  1.3  Ro  upstream  of  the  obsta¬ 
cle  nose  for  small  G  to  more  than  2.0  R*  for  the  largest 
G,  where  J2«  is  the  cometary  obstacle  nose  radius.  These 
trends  will  continue  as  still  larger  values  for  mass  loading 
are  considered,  or  similarly  if  smaller  obstacle  sizes  are 
assumed.  The  latter  possibility  highlights  a  significant 
remaining  need  for  a  model  for  the  flow  interior  to  the 
obstacle  boundary  to  locate  as  part  of  the  solution  the 
cometary  tangential  discontinuity  surface  that  separates 
the  flow  of  cometary  material  from  the  external  flow  of 
the  solar  wind  and  picked-up  ions. 

The  contour  lines  downstream  of  the  bow  wave  also 
change  dramatically  in  appearance  from  that  character¬ 
istic  of  planetary  flows  with  no  ion  pickup  at  the  smallest 
G  for  which  results  are  presented  to  being  nearly  concen¬ 
tric  semi-circles  (or  hemispheres)  about  the  obstacle  nose 
at  the  largest  G.  These  results  also  show  an  increasingly 
thick  dark  band  of  contour  lines  immediately  adjacent 


to  the  obstacle  indicating  an  enormous  enhancement  of 
the  density  in  that  region.  For  the  highest  mass  load¬ 
ing  determined  thus  far  in  this  study,  corresponding  to 
Sm  —  9.2,  our  results  indicate  densities  greater  than  40 
times  that  in  the  incident  solar  wind  at  the  largest  value 
for  G.  It  is  evident  there  exists  a  strong  tendency  for 
the  newly  ionised  particles  to  be  swept  back  toward  the 
obstacle  surface  where  they  form  a  dense  boundary  layer. 

Application  to  Venus 

Application  of  the  model  to  simulate  the  solar  wind 
interaction  with  Venus  requires  the  appropriate  source 
model.  For  this,  we  follow  the  suggestion  of  Nagy  et  al.", 
Breus",  and  Belotserkovskii  et  al."  and  assume  that  (a) 
the  only  important  source  of  new  ions  is  the  photoioniza¬ 
tion  of  the  hot  oxygen  corona,  (b)  any  initial  motion  of 
the  newly-created  ions  is  small  enough  that  there  is  no 
significant  addition  of  momentum  or  energy  to  the  flow, 
and  (c)  that  the  ion  source  term  can  be  represented  by 
an  expression  of  the  form 

Sm  =  Const  -  exp'  ^ T  (13) 

u„* 

in  which  J2o  is  the  distance  from  the  origin  at  the  plane¬ 
tary  center  to  the  nose  of  the  ionopause  obstacle,  and 
Conti  is  proportional  to  the  product  of  the  number 
density  no*  of  the  neutral  hot  oxygen  particles  at  the 
ionopause  nose  and  their  mass  m0* ,  divided  by  the  char¬ 
acteristic  time  r  for  photoionization.  Although  the  mass 
loading  beyond  the  bow  wave  may  be  small,  depending  on 
the  values  selected  for  H0+  and  Const  in  Equation  (13), 
in  the  results  reported  here,  we  have  not  terminated  the 
mass  addition  at  the  location  of  the  bow  wave  as  was 
done  in  Refs.  13  and  14,  but  allowed  its  influence  to  con¬ 
tinue  all  the  way  to  the  outer  boundary  of  the  grid.  In 
this  way,  our  solution,  including  in  particular  the  loca¬ 
tion  and  local  strength  of  the  bow  wave,  is  influenced  by 
essentially  all  of  the  mass  addition  beyond  the  ionopause 
and  not  just  that  downstream  of  the  bow  wave. 

Figure  6  presents  a  comparison  of  model  results  for  con¬ 
stant  density  contours  in  the  Venusian  ionosheath  with 
and  without  mass  loading  determined  as  described  above. 
The  calculations  were  carried  out  for  flow  conditions  cor¬ 
responding  to  those  of  orbit  582  of  the  Pioneer  Venus 
Or  biter  as  suggested  in  previous  studies  by  Mihalov  et 
al.",  Breus  ",  and  Belotserkovskii  et  al.14.  They  are 
that  the  oncoming  flow  has  a  free  steam  Mach  number 
of  5.7  and  a  value  of  5/3  for  the  ratio  of  specific  heats 
7,  and  the  Venusian  ionospheric  obstacle  has  a  shape  de¬ 
fined  by  a  value  of  0.03  for  the  ratio  h,/Rs  of  the  scale 
height  in  the  outer  ionosphere  to  the  obstacle  nose  dis¬ 
tance  from  the  planetary  center  in  accordance  with  the 
procedures4''’*  for  gasdynamic  modeling  of  the  flow  past 
Venus,  and  that  Const  and  H9+/R«  in  equation  (16)  are 
equal  to  4.49  and  0.062  respectively,  corresponding  to  a 
value  for  no*  of  3xl04  particles/cm',  H„+  =  400km  and 
Rs  =  6450km. 

Figure  6  shows  a  small  but  notable  difference  in  the 
density  contours  between  the  results  with  and  without 
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hot  oxygen  bum  loading  u  predicted  by  the  model.  The 
location  at  the  bow  wave,  however,  is  virtually  the  eame 
for  the  two  case*.  We  note  that  the  single  contour  line 
upetream  of  the  bow  wave  indicate*  tome  compression 
of  the  oncoming  solar  wind  to  occur  in  the  mass  loaded 
case.  The  significantly  sharp  bending  and  compaction 
of  the  density  contour  lines  near  the  ionopause  indicates 
the  presence  of  a  thin  boundary  layer  of  enhanced  density 
adjacent  to  the  ionopause. 

Figure  7  provides  a  comparison  of  the  result*  from  the 
current  model  and  that  of  Refs.  13  and  14  for  the  non¬ 
mass  loaded  and  mass  loaded  cases  for  Venus,  and  shows 
that  the  results  differ  substantantially  for  the  same  con¬ 
ditions.  Breus1*  and  Belotserkovski  et  al.14  find  the  mass 
loading  to  have  a  significant  effect  on  the  location  of  the 
bow  wave,  moving  it  to  a  position  slightly  beyond  that 
reported  by  Mihalov  et  al.1*,  whereas  our  results  indi¬ 
cate  almost  no  effect.  Their  results  also  indicate  the  bow 
wave  for  the  nonloaded  case  to  be  slightly  further  from 
the  ionopause  than  indicated  by  our  calculations,  with 
that  provided  in  the  more  recent  publication  of  Belot¬ 
serkovski  et  al.14,  presumably  corrected  from  those  re¬ 
ported  by  Breus1*,  being  in  better  agreement  with  our 
results.  The  shock  location  for  the  mass  loaded  case  as 
reported  in  both  Refs.  13  and  14  is  the  same.  Although 
that  result  is  in  better  agreement  with  .ae  observations 
of  Mihalov  et  al.1#,  we  have  carefully  reviewed  our  anal¬ 
ysis  and  believe  that  the  results  predicted  by  the  present 
model  are  correct  for  the  suggested  mas*  source  model 
and  mass  loading  parameters  suggested  for  this  example. 

To  determine  what  amount  of  mass  loading  would  be 
required  with  our  model  for  the  predicted  terminator  lo¬ 
cation  of  the  bow  wave  to  match  that  observed,  we  have 
carried  out  a  parametric  parameter  study  involving  an 
enhancement  of  the  mass  source  constant  (Const)  and 
the  hot  oxygen  scale  height  ( B a+)  as  given  in  the  source 
model  in  Equation  (16).  Figure  8  shows  the  results  for 
8  different  combinations  of  the  two  parameters.  In  the 
plot  on  the  left,  the  four  results  shown  for  the  bow  shock 
correspond  to  the  same  value  for  Const  as  used  for  Fig¬ 
ure  7,  namely  4.49  corresponding  to  a  number  density  of 
3  x  104  for  the  neutrals  at  the  ionopause  nose,  but  a  range 
of  values  from  400km  to  3200km  for  the  scale  height  B,+ 
at  the  hot  oxygen.  The  results  for  the  largest  of  these  val¬ 
ues  places  the  bow  wave  in  a  position  that  agrees  almost 
perfectly  with  the  observations  of  Mihalov  et  al.  (Ref. 
10),  but  the  value  of  3200km  seems  unreasonably  high 
for  the  scale  height  of  the  hot  oxygen. 

The  results  in  the  plot  in  the  right  half  of  Figure  8 
show  the  corresponding  results  for  four  other  combina¬ 
tions  of  the  two  mass  loading  parameters  (Const,  Bm*). 
They  show  that  the  bow  wave  moves  significantly  out¬ 
ward  from  the  planet  as  either  parameter  is  increased. 
Compared  with  these  four  bow  wave  locations,  the  data 
of  Mihalov  et  al.1*  for  Pioneer  Venus  Orbiter  for  orbit  S82 
falls  between  the  two  outermost  locations  corresponding 
to  na+  =  3  x  10*,  and  B,+  =  800  and  1600km.  It  appear* 
that*4  values  for  na*  of  this  order  may  be  appropriate  for 
Venus,  but  that  values  much  larger  than  about  400km  for 


the  scale  height  Bt+  at  the  hot  oxygen  are  probably  not 
appropriate.  A  substantial  part  of  the  outward  displace¬ 
ment  of  the  bow  wave  from  that  indicated  for  flows  with¬ 
out  mass  loading  thus  appears  to  be  accounted  for,  but 
some  additional  considerations  are  necessary  to  achieve 
a  good  quantitative  agreement  with  this  particular  data 
set.  Whether  the  main  source  of  these  differences  arises 
from  shortcomings  in  the  model  itself  or  in  the  interpre¬ 
tation  of  data  required  as  input  to  the  model  remains 
uncertain,  and  is  currently  under  study. 

All  of  the  numerical  calculations  for  which  results  are 
reported  here  were  performed  on  the  NASA  Ames  Re¬ 
search  Center  CRAY  XMP  computer  facility.  Approxi¬ 
mately  10  minutes  of  CPU  time  were  typically  required 
for  the  calculations  for  a  angle  case  of  the  more  difficult 
cometary  examples  considered  in  this  study,  and  approx¬ 
imately  2  minutes  for  the  typical  Venusian  mass  loaded 
cases.  Nearly  all  of  this  time  is  spent  with  the  nose  region 
solver,  as  the  tail  region  solver  is  extremely  rapid  and  can 
carry  the  solution  downstream  10  obstacle  nose  radii  in 
typically  less  than  30  CPU  seconds. 

Conclusions 

The  development  and  initial  verification  of  a  compu¬ 
tational  model  for  determining  the  solution  of  the  gas- 
dynamic  portion  of  the  gasdynamie  convected  magnetic 
field  model  for  solar  wind  flow  past  axisymmetric  obsta¬ 
cles  surrounded  by  an  extended  atmosphere  of  neutral 
particles  undergoing  ionization  so  as  to  add  mass,  mo¬ 
mentum,  and  energy  to  the  flow  is  described.  Results 
from  the  new  model  converge  continuously  to  those  cal¬ 
culated  previously  for  magnetoionospheres  without  ion 
pickup  using  completely  different  computational  proce¬ 
dures,  thereby  verifying  the  accuracy  of  both  methods. 

Specific  applications  are  made  to  cometary  interactions 
and  to  Venus  using  ion  loading  rates  provided  by  existing 
analyses  of  direct  observations.  For  comets,  the  high  mass 
loading  causes  the  bow  wave  to  be  weakened  due  to  the 
alowing  of  the  upstream  solar  wind  and  to  move  to  sub¬ 
stantially  greater  distances  from  the  nucleus.  For  Venus, 
the  ion  addition  is  comp&rtively  much  smaller,  and  the 
resulting  effects  more  modest,  although  still  significant 
for  quantitative  comparisons  between  theory  and  obser¬ 
vations.  In  both  applications,  large  density  enhancements 
are  indicated  by  the  model  immediately  adjacent  to  the 
obstacle  surface,  with  values  for  the  density  as  great  as 
40  times  that  in  the  incident  solar  wind  in  some  of  the 
cometary  examples.  Even  larger  enhancements  are  antic¬ 
ipated  when  the  model  is  applied  to  larger  comets,  such 
as  Halley,  with  more  massive  ion  additions  than  those  we 
have  considered  thus  far. 
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Figure  2.  Illustration  of  boundaries  and  coupling  of  com 
putational  domain*  in  the  mat*  loaded  solar  wind  plane 
tary/cometary  interaction  model. 


Figure  1.  Illustration  of  flow  field  subdivision,  extended 
grid,  and  coupled  flow  solver  combination  employed  in  the 
mass  loaded  solar  wind  planetary /cometary  interaction 
model. 
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Figure  3.  Illustration  of  a  typical  computational  grid  em¬ 
ployed  in  the  mass  loaded  solar  wind  planetary /cometary 
interaction  model. 


Figure  4.  Density  contours  for  a  simulated  cometary  flow 
with  extended  mass  loading  of  G  =  6  x  10”  particles /sec 
emitted  by  the  comet,  MSac  ~  6.0,  -t  =  5/3. 


Figure  6.  Calculated  d entity  contour*  in  the  Ve-mtian 
ionotheath  with  and  withour  matt  loading  due  to  pho¬ 
toionization  of  hot  oxygen,  Ms<x  *=  5.7,  7  *  5/3. 
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Figure  7.  Comparison  of  present  model  results  for  the 
location  of  the  Venusian  bow  shock  wave  with  those  of 
Breus13  and  Belotserkovski14  with  and  without  mass 
loading  due  to  photoionization  of  hot  oxygen,  Msoc,  — 
5.7,  7  =  5/3. 


Figure  8.  Predicted  locations  of  the  Venusian  bow  shock 
as  a  function  of  the  mass  source  model  constants  for  pho¬ 
toionisation  of  hot  oxygen,  M$ »  =  5.7,  7  =  5/3. 
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The  development  of  a  new  computational  model  for  predicting 
the  global  interaction  of  the  solar  wind  with  planetary 
magnetoionospheres  at  the  full  MHD  accuracy  level  for 
situations  when  the  IMF  is  aligned  with  the  oncoming  flow 
direction  is  described.  The  procedures  developed  are 
appropriate  for  supersonic,  super-Al fveni c  solar  wind  flows 
past  general  axisymmetric  magnetoionosphere  obstacle  shapes. 
The  computational  model  developed  employs  two  separate  but 
coupled  flow  field  solvers  to  determine  the  steady  flow  field: 
(1)  an  asymptotic  time-marching  procedure  for  the  subsolar  to 
terminator  region,  and  (2)  a  spatially-marching  procedure  for 
the  post-terminator  region.  Both  computational  procedures  are 
based  on  current,  fully-impl icit  numerical  algorithms,  and 
both  employ  fitted  discontinuity  representations  for  the  bow 
shock  and  magnetoionospheric  obstacle  surfaces.  An  outline  of 
the  numerical  model  and  some  of  the  computational  details 
involved  in  its  development  are  provided.  Results  of  initial 
applications  of  the  model  to  Earth  and  Venus  are  described. 
For  the  application  to  Earth,  the  results  are  compared  to 
those  previously  determined  by  Spreiter  and  Rizzi  (Acta  Astr., 
1974)  using  a  completely  different  computational  method.  The 
current  method  predicts  the  identical  variation  of  the  bow 
shock  shape  as  a  function  of  Alfven  Mach  number  in  the 
oncoming  solar  wind  that  was  previously  obtained,  i.e.  as  the 
oncoming  Alfven  Mach  number  decreases,  the  bow  shock 
simultaneously  moves  inward  toward  the  magnetopause  in  the 
nose  region  and  flares  outward  away  from  the  magnetopause 
along  its  flanks.  For  the  application  to  Venus,  the  model 
results  provide  an  important  indication  of  the  sensitivity  of 
the  Venusian  bow  shock  to  MHD  effects,  in  particular  the 
dependence  of  the  bow  shock  location  at  the  terminator  on  the 
Alfven  Mach  number  of  the  oncoming  flow. 
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Abstract 


The  gasdynamic  convected  magnetic  field  model  for  predicting  solar  wind 
flow  past  a  planetary  magnetoionopause  obstacle  has  been  extended  to 
three-dimensions  to  apply  to  obstacles  of  nonaxisymmetric  shape.  The 
need  for  such  an  extension  is  of  first-order  importance  for  Jupiter  and 
Saturn  because  the  effects  of  rapid  spin,  large  size,  and  substantial 
ring  current  phenomena  are  believed  to  cause  the  magnetospheres  of  these 
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planets  to  be  significantly  broader  near  the  planetary  equatorial  plane 
than  near  the  noon-midnight  polar  meridian  plane.  Direct  observational 
determination  of  such  asymmetries  for  Jupiter  and  Saturn  has  not  been 
possible,  however,  because  only  a  limited  amount  of  data  are  available 
from  four  spacecraft  at  Jupiter  and  three  at  Saturn,  all  of  which 
followed  flyby  trajectories  at  fairly  low  latitudes.  With  the  aid  of 
the  new  three-dimensional  model,  it  is  now  possible  to  infer  the  degree 
of  flattening  of  the  magnetospheric  cross  sections  from  a  knowledge  of 
the  locations  of  the  low-latitude  magnetopause  and  bow  shock  crossings. 
In  this  paper,  the  computational  procedures  of  the  new  model  are 
described,  and  calculated  results  are  presented  for  a  number  of  magneto¬ 
spheres  of  elliptic  cross  sections  with  values  ranging  from  1  to  2  for 
the  ratio  a/b  of  major  (equatorial)  to  minor  (polar)  axes.  This  range 
is  sufficient  to  include  values  appropriate  to  both  Jupiter  and  Saturn. 
Comparisons  of  the  model  results  with  the  locations  of  observed  cross¬ 
ings  of  the  magnetopause  and  bow  shock  directly  provide  an  estimate  of 
the  degree  of  equatorial  broadening  of  the  magnetospheric  cross  section 
for  each  of  these  planets.  For  Jupiter,  the  results  indicate  a 
broadening  to  a/b  ~  1.75,  a  value  consistent  with  previous  estimates 
determined  from  independent  calculations  of  the  three-dimensional 
magnetosphere  shape  formed  by  adding  to  the  planetary  dipole  field  the 
magnetic  field  of  an  equatorial  current  sheet  selected  so  as  to  match 
the  observed  Jovian  magnetic  field  near  the  equatorial  plane.  For 
Saturn,  a  similar  comparison  Indicates  a  smaller  broadening  to  a/b  ~ 
1.25.  The  determination  is  less  certain  than  for  Jupiter,  however, 
because  of  the  smaller  amount  and  greater  scatter  of  the  available  data. 
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Introduction 


The  gasdynamic  convected  magnetic  field  approximation  to  the  magne¬ 
tohydrodynamic  model  for  supersonic  solar  wind  flow  past  a  planetary 
magnetoionosphere  has  provided  a  useful  basis  for  understanding  the 
solar  wind  planetary  interaction  process  from  its  inception  more  than 
two  decades  ago  (Spreiter  et  al . .  1966)  to  the  present  day  (see  e.g., 
Spreiter  and  Stahara.  1985:  Slavin  et  al..  1985;  and  Russell.  1985,  for 
recent  reviews  of  the  theory  and  applications).  The  original  model  cal¬ 
culations.  although  accurate,  were  tedious  to  perform  with  the  methods 
available  at  the  time.  As  a  result,  solutions  were  carried  out  for  only 
a  limited  number  of  cases. 

More  recently  Spreiter  and  Stahara  (1980a, b)  completely  revised  the 
solution  procedures  without  changing  the  underlying  physics  in  order  to 
take  advantage  of  the  enormous  advances  in  computer  capability  and  com¬ 
putational  methods  that  have  occurred  since  the  original  model  was 
conceived.  They  produced  a  well -documented  model  (Stahara  et  al.,  1980) 
capable  of  routinely  providing  solutions  for  the  detailed  flow  and 
magnetic  field  properties  throughout  the  magnetoionosheath  region 
between  the  bow  shock  and  a  rather  arbitrarily  specified  smooth  shape 
for  the  planetary  magnetopause  or  ionopause.  The  results  provided  by 
the  model  are  now  widely  used  in  the  interpretation  of  spacecraft  data, 
and  In  further  extensions  of  the  theory. 

Throughout  the  theoretical  development  of  the  model,  the  shape  of  the 
magnetospheric  obstacle  has  always  been  approximated  as  axisymmetric 
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about  a  line  parallel  to  the  incident  aberrated  solar  wind  direction  and 

passing  through  the  planetary  center.  For  the  earth,  this  has  been 

justified  by  numerous  studies  of  extensive  observational  data  that  have 

shown  nat  any  departure  from  axisymmetry  is  virtually  immeasurable. 

More  precisely.  Fairfield  (1971)  and  Holzer  and  Slavin  (1978)  have 

determined  that  the  eccentricity  e  of  the  magnetospheric  cross  section 

in  the  terminator  plane  Is  less  than  0.2.  In  this,  the  eccentricity  Is 

2  2  1/2 

defined  in  the  usual  way  as  the  positive  root  of  t  -  (l*b  /a  )  where 

a  and  b  are  the  major  and  minor  axes  of  an  ellipse  that  best  fits  the 

observed  magnetospheric  cross  section  in  the  terminator  plane.  Since  e 

<  0.2  correspond  to  values  for  a/b  between  1  and  1.02  as  can  be  readily 

2  *1/2 

determined  using  the  equivalent  relation  a/b  -  (1-  e  )  ,  the  approxi¬ 

mating  elliptic  cross  section  for  the  earth’s  magnetosphere  differs  only 
slightly  from  being  circular.  Furthermore,  predicted  results  based  upon 
an  axisymmetric  magnetopause  shape  have  consistently  been  found  to  be  in 
good  agreement  with  observations  of  both  the  bow  shock  location  and  also 
of  the  plasma  and  magnetic  field  properties  throughout  the  magnetosheath 
of  the  earth.  Similar  results  have  been  found  for  the  other  terrestrial 
planets  (Slavin  et  al..  1983),  but  notable  differences  appear  in  the 
corresponding  comparisons  for  Jupiter  and  Saturn  (Slavin  et  al..  1985). 

All  direct  measurements  of  conditions  In  the  magnetospheres  and 
surrounding  solar  wind  Interaction  regions  at  Jupiter  and  Saturn  are 
those  provided  by  flyby  missions  of  Pioneer  and  Voyager  spacecraft. 
Their  trajectories  as  well  as  the  average  observed  locations  of  the  bow 
shock  and  magnetopause  boundaries  are  shown  in  Figure  1  (see  Slavin  et 
al..  (1985)  and  references  therein).  The  spatial  coverage  provided  by 
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the  four  spacecraft  trajectories  through  the  Jovian  system  and  the  three 
through  the  Saturnian  system  appears  to  be  reasonably  good,  as  viewed  in 
the  rotated  cylindrical  coordinate  system  of  Figure  1.  However,  all  of 
the  trajectories  were  at  fairly  low  latitudes,  and  no  data  whatsoever 
were  obtained  at  high  latitudes. 

In  spite  of  the  scatter  in  the  average  boundary  and  bow  shock  locations, 

much  of  which  can  be  accounted  for  by  the  effects  of  differing  solar 
wind  stagnation  pressure  p  ,*  p  U  ,  an  orderly  pattern  of  boundary 

and  shock  crossings  at  Jupiter  can  readily  be  recognized  in  Figure  1 
even  before  application  of  corrections  for  variations  in  p  .  For 

w  t 

Saturn,  the  situation  is  not  so  clear  because  the  data  display  more 

scatter  than  can  be  accounted  for  by  a  simple  scaling  of  the  magneto¬ 
sphere  and  bow  shock  according  to  known  variations  in  p$t..  It  has  been 

suggested,  furthermore,  that  the  distant  locations  of  the  magnetopause 
and  bow  shock  observed  by  Voyager  2  may  be  the  result  of  Saturn  possibly 
being  in  the  distant  magnetic  tail  of  Jupiter  at  the  time  of  the 
encounter  (see  e.g..  Scarf  et  al.,  1977;  Grzedzielski  et  al.,  1981; 
Behannon  et  al . ,  19bJ;  and  Oesch,  1983  for  further  discussion  of  this 
interesting  possibility).  Most  of  the  expected  effects  of  such  an 
occurrence  should  be  compensated  for,  however,  by  the  stagnation  pres¬ 
sure  scaling  procedure  of  Slavln  et  al .  (1985)  that  has  been  applied  to 
the  observational  data  summarized  In  Figure  2,  and  subsequently  herein. 

In  a  comprehensive  study  of  the  positions  of  the  magnetopause  and  bow 
shock  of  Jupiter  and  Saturn.  Slavin  et  *1.  (1985)  determined  the  least- 
square  conic  curve  representations  recast  into  nondimensional  form  in 
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Figure  2.  The  results  for  Jupiter  show  that  the  calculated  bow  shock 
based  upon  an  assumed  axisymmetric  magnetopause  shape  that  matches  the 
observed  equatorial  magnetopause  Is  everywhere  displaced  substantially 
farther  from  the  magnetopause  than  the  observations.  The  results  for 
Saturn  display  a  similar  result  near  the  nose  of  the  bow  shock. 
However,  both  the  original  least-square  conic  curve  of  Slavin  et  al . 
(1985),  and  an  analogous  new  curve  that  includes  effects  of  outbound 
Voyager  1  data  excluded  in  the  original  determination  may  be  seen  to 
flare  out  even  beyond  the  calculated  axisymmetric  shock  near  or 
downstream  of  the  terminator  plane.  Temporarily  setting  aside  this 
flaring  at  Saturn,  which  will  be  discussed  separately,  the  differences 
between  the  calculated  axisymmetric  and  the  least-square  observational 
shock  locations  in  the  nose  region  of  each  of  these  planets  are 
substantial.  More  specifically,  the  calculated  thickness  of  the 
magnetosheath  in  the  subsolar  region  is  80  percent  too  great  for  Jupiter 
and  30  percent  too  great  for  Saturn. 

Slavin  et  al .  (1985).  attributed  these  differences  to  a  significant 
departure  from  axisymmetry  in  the  magnetospheri c  shapes  for  these 
planets.  In  this  view,  the  equatorial  regions  of  the  magnetopause  are 
broadened  with  respect  to  the  polar  regions  by  effects  of  a  large  equa¬ 
torial  ring  current  system  (Engle  and  Beard,  1980)  brought  about, 

presumably,  by  centrifugal  effects  associated  with  the  large  angular 
velocity  <o  and  low  mean  density  of  these  planets. 

The  equatorial  broadening  hypothesis  Is  consistent  with  most  of  the 
results  shown  in  Figure  2  since  the  axisymmetric  shape  that  results  from 
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rotating  the  observed  equatorial  magnetopause  trace  about  the  symmetry 
axis  produces  a  shape  that  has  a  significantly  larger  cross-section  area 
than  a  flattened  shape  having  the  same  equatorial  magnetopause  trace. 
As  a  result,  the  bow  shock  about  an  axisymmetric  magnetopause  will  be 
everywhere  farther  from  the  magnetopause  than  it  Is  for  a  flattened  mag¬ 
netopause  having  the  same  equatorial  shape,  just  as  shown  for  Jupiter 
and  near  the  magnetopause  nose  for  Saturn. 

To  appreciate  the  enhanced  effects  of  centrifugal  forces  in  the  Jupiter 

and  Saturn  magnetospheres  compared  with  those  for  the  earth,  note  that 
the  synchronous  orbital  distance  r$  beyond  which  the  centrifugal  force 
2 

rco  per  unit  mass  corotating  with  the  planet  in  the  equatorial  plane 

2  3  2  -fi 

exceeds  the  gravitational  force  GM/r  -4jtGppRp  /3r  ,  where  G-6. 668x10 
2  -  2 

dynes  cm  g  is  the  gravitational  constant,  and  Rp  is  the  planetary 

2  1/3 

radius,  is  at  a  normalized  planetary  distance  rs / Rp  -  [4jiGp/(3o)  )]  ~ 

6.6  for  the  earth,  but  only  about  2.2  and  1.8  for  Jupiter  and  Saturn. 

However,  a  more  appropriate  indicator  of  centrifugal  flattening  of  the 
magnetosphere  is  the  ratio  of  rg  to  a  representative  dimension  of  the 

magnetospheric  obstacle,  rather  than  the  planetary  radius.  If  this 
dimension  is  taken  to  be  the  planetocentric  distance  RQ  to  the  magne¬ 
topause  nose,  where  Rp  -  (10Re,  58.6Rj.  19. 7R^)  for  the  earth,  Jupiter 
and  Saturn,  the  resulting  values  for  r  /R  are  0.66  for  the  earth,  and 

5  0 

only  0.04  and  0.09  for  Jupiter  and  Saturn.  This  clearly  shows  that 
centrifugal  force  on  corotating  matter  dominates  over  gravitational 
force  throughout  a  much  larger  fraction  of  the  Jupiter  and  Saturn  magne¬ 
tospheres  than  it  does  for  the  earth. 
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On  the  basis  of  the  available  observations,  Slavin  et  al.  (1985) 
proposed  an  equatorial  broadening  hierarchy  similar  to  that  illustrated 
in  Figure  3  for  the  magnetopause  cross  sections  in  the  terminator  plane. 
As  discussed  above,  the  cross  section  of  the  earth's  magnetosphere  is 
essentially  circular  with  an  eccentricity  0.2,  which  corresponds  to 
a  ratio  of  major  to  minor  axes  of  1.00  £  a/b  £  1.02.  In  contrast,  they 
concluded  that  Jupiter's  magnetosphere  Is  broadened  Into  an  ellipse  with 
e  -  0.8,  or  a/b  -  5/3.  The  corresponding  values  for  Saturn  were  not 
determined  with  an  equal  precision,  but  it  was  concluded  that  planet 
presents  an  intermediate  case  with  0  £  e  £  0.8,  or  1  <  a/b  £  5/3.  Our 
results  are  most  consistent  with  the  Pioneer  and  Voyager  observations 
when  values  for  a/b  are  set  to  1.25  for  Saturn  and  1.75  for  Jupiter,  and 
the  cross  sections  for  those  planets  are  drawn  accordingly  in  Figure  3. 
In  each  case,  the  shock  position  is  that  calculated  at  the  terminator 
plane  using  representative  values  of  2  for  the  ratio  of  specific  heats 
Y  ,  and  8.  10  and  12  for  the  Mach  number  of  the  solar  wind  approaching 
earth.  Jupiter  and  Saturn. 

To  provide  the  appropriate  theoretical  solutions  for  application  to 
equatorially  broadened  magnetospheres,  we  have  extended  our  previous 
axisymmetric  analysis  so  that  it  can  be  applied  to  the  interaction  of 
the  solar  wind  flow  with  a  general  three-dimensional  magnetoionospheric 
obstacle.  This  extension  serves  the  dual  purpose  of  providing  (a)  more 
accurate  theoretical  capability  required  for  comparative  studies  with 
observations  of  the  solar  wind  interaction  with  substantially  nonaxisym- 
metric  planetary  magnetospheres,  and  (b)  an  important  element  of  a 
general  interactive  scheme  in  which  the  three-dimensional  shape  of  the 
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magnetoi onopause  and  the  surrounding  flow  and  magnetic  fields  for  any 
planet  are  calculated  alternately  in  the  process  of  converging  toward  an 
exact  solution  of  the  magnetohydrodynamic  equations. 

The  new  computational  model  is  described  in  this  paper,  and  results  are 
presented  for  the  calculated  shape  of  the  bow  shock  and  thr.  flow 
properties  in  the  magnetosheaths  of  equatorially  broadened  magneto¬ 
spheres  having  elliptic  cross  sections  representative  of  those  suggested 
for  Jupiter  and  Saturn.  By  comparison  of  these  results  with  the 
observed  locations  of  the  low-latitude  magnetopause  and  bow  shock  cross¬ 
ings.  a  quantitative  determination  is  obtained  for  the  degree  of  equato¬ 
rial  broadening  of  the  magnetospheres  of  Jupiter  and  Saturn.  The  new 
three-dimensional  results  for  these  broadened  magnetospheres  are  shown 
to  be  in  satisfactory  agreement  with  the  observations  near  these 
planets,  and  to  account  for  most  of  the  differences  between  the  observa¬ 
tions  and  the  axisymmetric  model  results  indicated  previously  by  Slavin 
et  al.  (1985). 


Mathematical  Model  and  Solution  Procedure 

The  fundamental  assumption  underlying  the  present  work,  as  in  our 
previous  analyses  of  axisymmetric  flows  reported  In  the  references  cited 
above,  is  that  the  bulk  properties  of  the  solar  wind  plasma  flow  around 
a  planetary  magnetoionosphere  can  be  described  adequately  for  many  pur¬ 
poses  by  solutions  of  the  continuum  equations  of  magnetohydrodynamics  of 
a  perfect  gas  having  infinite  electrical  conductivity  and  zero  viscosity 
and  thermal  conductivity.  The  primary  justification  for  this  assumption 
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is  provided  not  by  rigorous  proof,  but  by  the  outstanding  general  agree¬ 
ment  of  the  calculated  results  with  a  large  amount  of  direct  observa¬ 
tions  in  space. 

Under  this  assumption,  the  equations  used  to  describe  the  flow  of  the 
solar  wind  plasma  past  a  planetary  magnetoionosphere  are  those  of 
single-fluid  dissipationless  magnetohydrodynamics  for  the  conservation 
of  mass,  momentum,  energy,  and  magnetic  flux,  augmented  by  the  equation 
of  state  for  a  perfect  gas  and  by  entropy  considerations.  These  equa¬ 
tions  have  been  given  many  times  (see  Spreiter  and  Stahara.  1985,  for 
our  most  recent  account),  and  will  not  be  repeated  here. 

In  addition  to  the  differential  equations  and  auxiliary  relations 
referred  to  above,  boundary  and  initial  conditions  must  be  specified  for 
the  properties  of  the  solar  wind  upstream  of  the  bow  shock  wave,  and  for 
the  size  and  other  relevant  properties  of  the  planet  for  which  a  solu¬ 
tion  is  sought.  For  a  planet  with  a  strong  dipolar  component  to  the 
magnetic  field,  like  the  earth,  Jupiter,  and  Saturn,  the  pertinent 
parameters  are  the  strength  and  orientation  of  the  planetary  field. 

In  common  with  most  studies  of  the  present  type,  we  assume  that  the  flow 
Is  steady  and  that  the  dipole  field  is  stationary  in  time.  As  will  be 
commented  upon  further  in  the  discussion  of  the  results,  the  large 
dimensions  of  the  magnetospheres  of  Jupiter  and  Saturn  cause  these 
assumptions  to  be  much  more  restrictive  than  they  are  in  applications  to 
the  earth  and  other  terrestrial  planets. 


€4 


For  typical  solar  wind  conditions,  particularly  at  the  locations  of  the 
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outer  planets,  both  the  oncoming  sonic  Mach  number  Mc  -  V  / ( y p  /p  ) 
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and  the  Alfven  Mach  number  M.  -  V  / ( B  */(4np  l)1^  are  high  (of  the 

order  of  10).  Under  these  conditions,  an  important  simplification  of 
the  magnetohydrodynamic  equations  occurs  because  the  magnetic  terms  may 
be  disregarded  as  small  in  the  momentum  and  energy  conservation 
equations.  With  these  terms  omitted,  the  magnetohydrodynamic  equations 
separate  into  two  sets.  One  is  the  ordinary  gasdynamic  equations  for 
the  flow  of  a  perfect  inviscid  gas.  The  other  set  is  the  Faraday 
magnetism  equations  for  a  perfectly  conducting  medium. 

The  equations  for  the  fluid  motion  about  an  obstacle  of  given  shape  thus 
reduce  to  those  of  gasdynamics.  and  can  be  solved  without  further 
consideration  of  the  magnetic  field.  The  magnetic  field  can  be  deter¬ 
mined  subsequently  by  solving  the  remaining  Faraday  equations  using  the 
values  for  the  velocity  and  density  fields  calculated  using  the  gasdy¬ 
namic  equations.  Because  the  magnetic  field  determined  in  this  fashion 
is  usually  interpreted,  although  somewhat  ambiguously,  as  being  "frozen- 
in"  or  "moving  with  the  fluid",  the  resulting  equations  and  the  model 
they  represent  are  commonly  referred  to  as  the  gasdynamic  convected 
field  approximation.  It  should  be  noted,  however,  that  the  interpreta¬ 
tion  of  the  field  as  frozen-in  or  moving  with  the  fluid  is  not  dependent 
upon  the  separation  of  the  determination  of  the  fluid  motion  from  that 
of  the  magnetic  field,  but  Is  actually  an  exact  consequence  of  the  basic 
assumption  that  the  flow  can  be  represented  by  the  magnetohydrodynamic 
equations  for  a  dissipationless  electrically  conducting  perfect  gas. 
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In  carrying  out  the  decoupled  gasdynamic  calculation,  the  shape  of  the 
magnetoionopause  is  not  known  a  priori,  but  must  be  found  as  part  of  the 
overall  solution.  In  our  previous  analyses,  the  magnetoionopause  shape 
has  been  determined  by  either  an  independent  approximate  theoretical 
analysis,  a  synthesis  of  actual  observational  results,  or  a  combination 
thereof.  The  latter  procedure  has  also  been  employed  in  the  present 
analysis. 

Numerical  solutions  of  the  full  magnetohydrodynamic  equation  set  without 
resort  to  the  gasdynamic  convected  field  approximation  now  exist,  but 
they  are  still  in  a  fairly  early  state  of  development,  require  large 
computational  resources,  and  lack  adequate  resolution  for  many  purposes. 
In  summarizing  the  situation  recently.  Russell  (1985)  stated. 
"Currently,  the  only  way  to  achieve  the  spatial  resolution  needed  for  a 
useful  comparison  with  data  obtained  on  a  pass  through  a  planetary 
magnetosheath  is  to  employ  the  gasdynamic  convected  field  model". 

One  of  the  most  important  features  of  both  the  current  three-dimensional 
and  previous  axisymmetric  models  is  the  use  of  two  separate,  but 
coupled,  computational  procedures  for  determining  the  steady  flow  field. 
As  illustrated  In  Figure  4,  the  first  of  these  is  a  nose  region  solver 
that  determines  the  transonic  flow  field  from  the  subsolar  region  to  or 
somewhat  beyond  the  terminator  location  by  using  an  unsteady  procedure 
that  marches  the  solution  forward  in  time  from  a  somewhat  arbitrarily 
imposed  initial  state  to  the  asymptotic  steady  state  solution.  That 
solver  is  then  linked  together  with  a  spatially-marching  tail  region 
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solver  that  advances  the  supersonic  flow  field  solution  downstream  to 
any  arbitrary  distance  in  which  the  flow  remains  supersonic. 

The  choice  of  axial  location  at  which  the  two  solutions  are  joined  is 
relatively  arbitrary,  the  most  important  requirement  being  that  Imposed 
by  the  spatially-marching  solver  that  the  axial  component  of  the  local 
sonic  Mach  number  be  greater  than  unity  everywhere  at  the  joining  plane. 
In  most  of  our  previous  applications,  it  was  convenient  to  join  the  two 
solutions  at  the  plane  through  the  planetary  center  oriented  normal  to 
the  oncoming  solar  wind  direction,  which  we  term  the  terminator  plane. 
For  some  cases  involving  magnetoi onopause  obstacles  with  enhanced 
flaring  at  the  terminator,  it  has  been  found  necessary  to  join  the 
solutions  at  a  plane  some  fraction  of  an  obstacle  nose  radii  downstream 
of  the  terminator  in  order  to  satisfy  the  requirement  of  supersonic 
axial  velocity  component. 

For  clarification,  it  should  be  noted  that  our  definition  of  the 
terminator  plane  differs  slightly  from  the  standard  definition  based  on 
visual  considerations  as  the  plane  through  the  planetary  center  oriented 
normal  to  the  sun-planet  line.  The  orientations  of  the  two  planes 
differ  by  the  small  angle  between  the  upstream  direction  and  the  sun- 
planet  line,  and  is  due  primarily  to  the  aberration  effects  of  the 
orbital  motion  of  the  planet  through  the  much  higher  speed  and  nearly 
radial  outflow  of  the  solar  wind  from  the  Sun.  For  earth,  this  angle  is 
usually  about  5°.  but  substantially  smaller  values  would  be  anticipated 
for  Jupiter  and  Saturn.  The  reason  is  that  (a)  the  average  radial 
velocity  of  the  solar  wind  is  nearly  independent  of  distance  from  the 
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Sun  throughout  this  part  of  the  Solar  System,  and  (b)  the  orbital  velo¬ 
city  of  a  planet  varies  inversely  as  the  square  root  of  its  distance 
from  the  Sun.  Since  this  distance  Is  approximately  5.2  and  9.5  times 
greater  for  Jupiter  and  Saturn  than  for  the  earth,  the  corresponding 
difference  between  the  orientations  of  the  alternately  defined  termina¬ 
tor  planes  is  only  about  2.2°  for  Jupiter  and  1.6°  for  Saturn. 

The  reason  for  expending  the  substantial  additional  effort  to  develop 
such  a  coupled  two-flow  solver  combination,  rather  than  employing  a 
single  time-marching  solver  for  the  entire  flow  field  as  is  done  in  all 
existing  magnetohydrodynamic  solutions  procedures,  are  twofold.  First, 
there  is  a  huge  gain  in  computational  efficiency  because  the  com¬ 
putationally  expensive  time-marching  procedure  required  to  treat  the 
subsonic/transonic  flow  in  the  nose  region  is  used  only  where  needed, 
while  the  remainder  of  the  solution  is  obtained  using  a  highly  efficient 
spatially  marching  procedure  appropriate  for  supersonic  flow.  Secondly, 
a  computational  capability  is  attained  to  determine  flow  conditions  as 
arbitrarily  far  downstream  alongside  a  specified  magnetotail  surface  as 
needed.  Such  an  extended  capability  cannot  be  provided  by  a  time¬ 
marching  procedure  alone.  Finally,  we  note  that  such  a  subdivision  of 
the  solution  procedures  is  quite  common  in  calculating  super¬ 
sonic/hypersonic  flows  past  blunt  nosed  bodies.  It  has  been  used 
throughout  not  only  all  our  previous  work,  but  also,  for  example,  in  the 
analysis  of  flow  past  a  reentry  body  such  as  the  space  shuttle. 
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In  the  nose  region,  we  have  adapted  the  procedure  of  filzk  et  al.  (1963) 
to  solve  the  unsteady,  three-dimensional  dissipationless  gasdynamlc 
partial  differential  equations  governing  the  conservation  of  aass, 
■omentum,  and  energy.  Upon  Introduction  of  the  generalized  Independent 
variable  transformation  t  -  t.  $  -  &(x.y.z).  -  ti(x.y.z),  C  “ 
C(x.y.z).  these  equations  can  be  written  In  strong  conservation  law  form 
as: 

where 
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Here  (u.v.w)  are  the  Cartesian  velocity  components,  p  and  p  are  density 

and  pressure,  e  is  the  internal  energy  which  is  related  to  pressure  by  p 
2  2  2 

-  (y-l)[e-  (u  +  v  +  w  )/2],  where  y  is  the  ratio  of  specific  heats. 
(U.V.W)  are  the  contravariant  velocities  without  metric  normalization 
given  by 
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in  wnich  (£x>  \  Sz).  (tix.  t^.  ) .  (£x.  Cz>  are  the  transformation 

metrics,  and  J  is  the  Jacobian  of  the  transformation.  At  each  time  step 
the  transformation  metrics  are  unknown  and  must  be  determined  as  part  of 
the  solution.  In  our  implementation,  they  are  evaluated  numerically 
usirj  second-order  central  differences.  The  time  dependent  metric  terms 
are  functions  of  the  instantaneous  shock  speed,  which  is  determined  from 
the  shock  fitting  procedure  of  Kutler  et  al.  (1980). 

The  numerical  algorithm  used  to  solve  the  above  equations  is  based  on 
the  Beam  and  Warming  (1976)  procedure.  The  method  is  employed  in  delta 
form  and  is  fully  Implicit,  noniterative,  second-order  accurate  in 
space,  and  spatially  factored.  When  applied  to  equation  (1),  the  algo¬ 
rithm  can  be  written  as: 
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AAA  A  AA^  AAA  AA 

where  (A.B.C)  are  the  Jacobian  matrices  A  -  dE/Jq,  B  -  dF/dq.  C  -  dG/dq. 

A  fl+l  fl 

Aq  -  (q/J)  -  (q/J)  where  n  represents  the  time  level,  (8^.  8^.  8^) 

represent  second-order  central  difference  operators,  (e^.  eg)  are 

implicit  second-order  and  explicit  fourth-order  smoothing  parameters, 
( V, a)  represent  first-order  backward  and  forward  difference  operators, 
and  t  is  the  time  step.  A  single  time  step  integration  of  equation  (4) 
consists  of  three  separate  inversion  steps  associated  with  each  of  the 
three  spatial  directions  ( 4 . "H •  C > -  Each  inversion  step  involves  the 
solution  of  a  5x5  block  tridiagonal  system.  Integration  step  size  is 
established  by  using  the  maximum  eigenvalue  of  the  Jacobian  matrices 

AAA 

(A.B.C). 

The  analysis  is  initiated  by  introducing  a  computational  mesh  to 
discretize  the  flow  field.  For  the  present  three-dimensional  model  as 
well  as  in  all  our  previous  work,  we  employ  a  fitted  discontinuity  sur¬ 
face  representation  of  both  the  bow  shock  and  the  magnetospheric 
obstacle.  For  application  to  three-dimensional  planetary  magnetospheres 
we  have  found  It  convenient  to  develop  a  mesh  generator  based  upon  a 
generalized  spherical  (r,  0,  0)  coordinate  •system  with  the  origin 
usually  located  at  the  planetary  center,  but  occasionally  shifted 
downstream  from  that  location  on  a  line  through  the  planetary  center  and 
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parallel  to  the  oncoming  solar  wind.  Although  the  basic  method  imposes 
no  such  requirement  for  the  calculations  presented  here,  it  has  been 
assumed  that  the  magnetopause  shape  and  associated  flow  field  possess 
bilateral  or  quadrilateral  symmetry  about  the  equatorial  plane  and  the 
polar  meridian  plane. 

The  associated  Cartesian  (x,y,z)  coordinate  axes  and  a  typical  grid 
employed  in  our  analysis  are  illustrated  in  Figure  5  where  a  perspective 
cut  out  view  of  one  quadrant  of  the  flow  field  is  provided.  In  that 
figure  the  origin  coincides  with  a  location  approximately  half  of  an 
obstacle  radius  downstream  of  the  planetary  center  and  the  x*axis.  which 
coincides  with  the  longitudinal  axis  of  the  magnetoionosphere,  points 
directly  into  the  oncoming  solar  wind.  The  y  axis  lies  in  the  planetary 
magnetic  equatorial  plane,  and  the  2  axis  points  in  the  north  polar 
direction.  In  Figure  5.  the  obstacle  to  shock  (r.  6)  grids  in  the 
azimuthal  planes  containing  the  x-z  <0  -  0°)  and  x-y  (0  -  90°)  axes  are 
shown  together  with  the  mesh  distribution  on  the  magnetoi onopause 
surface,  and  the  outflow  boundary  grid  located  in  the  y*z  plane  at  the 
final  downstream  x  location. 

The  quadrant  of  the  (r,  e.  $)  mesh  shown  in  Figure  5  contains  (19,30,9) 
points,  with  the  9  azimuthal  planes  spaced  at  equal  10°  increments.  In 
each  azimuthal  plane,  the  capability  of  generating  a  body  normal  mesh 
has  been  Implemented  In  the  grid  generator.  A  clustering  capability  has 
also  been  Implemented  to  enable  Independent  grid  point  clustering  in 
each  of  the  three  coordinate  directions.  For  the  magnetospher i c 
obstacle  shown  in  Figure  5.  which  represents  an  equatorial  broadened 
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Jupiter  magnetopause  with  a/b  -  1.75,  the  equal  angular  distribution  of 
the  azimuthal  planes  results  in  a  natural  clustering  of  the  grid  spacing 
on  the  body  in  the  polar  region  In  contrast  to  that  In  the  equatorial 
region.  This  can  be  retained  or  easily  modified  If  desired  by  employing 
the  clustering  option. 

Boundary  conditions  necessary  for  specifying  a  properly  posed  math¬ 
ematical  problem  are  that  the  flow  satisfy  the  unsteady  Rankine-Hugoniot 
shock  relations  along  the  moving  bow  shock  surface,  be  entirely  super¬ 
sonic  at  the  downstream  outflow  boundary,  possess  bilateral  or  quadri¬ 
lateral  symmetry  about  the  stagnation  streamline  at  the  nose,  and  be 
tangential  to  the  obstacle  at  its  surface.  Initial  flow  field  condi¬ 
tions  are  determined  by  use  of  an  approximating  surface  for  the  bow 
shock  wave,  and  by  prescribing  a  modified  Newtonian  pressure  distribu¬ 
tion  on  the  surface  of  the  magnetoi onopause .  Since  the  maximum  entropy 
streamline  wets  the  obstacle,  that  fact  plus  the  flow  tangency  condition 
on  the  magnetoi onopause  serve  to  determine  the  remainder  of  the  flow 
properties  on  that  surface.  A  linear  variation  for  the  flow  properties 
between  the  bow  shock  and  the  magnetoionopause  then  provides  the  start¬ 
ing  flow  field  which  Is  then  Integrated  in  a  time-asymptotic  fashion 
until  the  steady  state  solution  is  obtained. 

At  the  boundaries,  modification  of  the  differencing  algorithm  to  account 
for  the  appropriate  conditions  described  above  Is  accomplished  as 
follows.  For  the  fitted  discontinuity  surface  representing  the  bow 
shock,  the  flow  variables  immediately  downstream  of  the  shock  are  deter¬ 
mined  using  the  unsteady  shock  fitting  procedure  of  Kutler  et  al. 
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(1980).  At  the  outflow  boundary  where  the  flow  is  entirely  supersonic, 
the  dependent  variables  are  determined  by  extrapolation  from  the  adja¬ 
cent  interior  points.  Along  the  stagnation  streamline,  which  forms  a 
singularity  line  due  to  our  use  of  a  spherical -1 i ke  coordinate  system 
(Pulliam  and  Steger,  1978).  symmetry  conditions  are  enforced  by  employ¬ 
ing  a  second-order-accurate  averaging  procedure.  Also  incorporated  in 
the  numerical  solver  is  the  capability  to  enforce  symmetry  boundary 
conditions  across  any  two  arbitrary  azimuthal  planes.  For  the  present 
application,  all  members  of  the  family  of  three-dimensional  obstacles 
being  considered  possess  quadrilateral  symmetry  about  two  orthogonal 
azimuthal  planes  passing  through  the  equator  and  the  polar  meridians. 
Consequently,  the  flow  calculations  can  be  restricted  to  a  single  quad¬ 
rant  as  shown  in  Figure  5.  thereby  reducing  by  a  factor  of  four  the 
number  of  grid  points  at  which  the  solution  must  be  determined.  The 
obstacle  surface  flow  tangency  condition  is  incorporated  by  setting  the 
contravari ant  velocity  component  normal  to  the  surface  equal  to  zero, 
and  then  determining  the  other  two  contravariant  velocity  components  and 
the  density  at  the  obstacle  surface  by  extrapolation  from  interior 
points.  The  total  energy  is  then  evaluated  by  using  that  information 
together  with  the  pressure  computed  from  the  normal  momentum  equation. 
The  interior  flow  field  bounded  by  these  various  boundaries  is  treated 
In  shock-capturing  fashion  and  therefore  allows  for  the  accurate  deter¬ 
mination  of  secondary  shocks,  should  any  occur. 
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In  the  tall  region,  we  have  employed  the  shock  capturing  technique  of 
Kutler  et  al.  (1973)  to  solve  the  steady,  three-dimensional  dlsslpa- 
tlonless  gasdynamlc  equations.  The  conservation  equations  for  mass  and 
momentum  written  In  weak  conservation  law  form  In  cylindrical  (x.R.$) 
coordinates  can  be  written  as: 

®x+^R+^+R*°  (5) 

where  the  four-component  vectors  <0,F,(a,H)  are  defined  by 


In  which  (u.v.w)  denote  velocity  components  In  the  (x,R,$)  coordinate 
directions.  The  governing  set  of  equations  Is  completed  with  the  addi¬ 
tion  of  energy  conservation  as  given  by  the  equation  for  the  total 
enthalpy 
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Ht  -  h(p,p)  +  q  /2  -  const 


Here  q  is  the  magnitude  of  the  velocity  vector  and  h(p.p)  represents  the 
static  enthalpy  which  for  a  perfect  gas  Is  given  by  [y/(y-l ) 3 ( p/p ) . 

The  governing  equations  are  transformed  to  a  computational  space  using 
the  independent  variable  transformation  %  -  x.  -  (R*Rb)/(Rj-Rb>.  C 

-  4  .  where  (R^.R^)  represent  the  cylindrical  radial  coordinate  of  the 

bow  shock  and  body  surface,  respectively,  and  $  Is  the  azimuthal  angle. 
This  yields  the  conservation  equation 
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The  numerical  algorithm  chosen  to  solve  equation  (8)  is  MacCormack’s 
explicit  second-order  accurate  predictor-corrector  method.  As  applied 
to  equation  (9),  the  algorithm  can  be  written  as: 
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Subsequent  to  each  Integration  step  of  equation  (10)  with  respect  to  the 

hyperbolic  marching  coordinate  the  physical  flow  variables 

(p.p.u.v.w)  must  be  decoded  from  the  components  ui  of  U.  This  necessi- 
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tates  the  solution  of  five  simultaneous  nonlinear  equations  consisting 

2 

of  equation  (7)  together  with  the  four  elements  ui  -  (pu,  p+pu  .  puv. 

puw).  For  a  perfect  gas.  the  procedure  yields  a  quadratic  equation  that 
can  be  solved  for  the  supersonic  velocity  component  u.  The  remaining 
physical  variables  can  then  be  determined  directly. 

At  the  bow  shock,  the  fitted  shock  discontinuity  approach  of  Thomas  et 
al .  (1972  )  is  employed,  while  at  the  obstacle  surface  the  surface 
tangency  condition  is  enforced  by  applying  the  method  of  Abbett  (1971). 
Integration  step  size  in  the  marching  direction  is  redetermined  at  each 
location  so  as  to  allow  the  maximum  possible  step  size  consistent  with 
stabi 1 ity. 

The  numerical  calculations  for  which  results  are  reported  here  were 
performed  on  the  NASA  Ames  Research  Center  CRAY  XMP  computer  facility. 
Approximately  20  to  30  minutes  of  CPU  time  were  typically  required  for 
the  calculations  for  a  single  case  of  the  three-dimensional  magnetoiono¬ 
sphere  shapes  considered  in  this  study.  Essentially  all  of  this  time  is 
spent  with  the  nose  region  solver.  For  the  family  of  three-dimensional 
shapes  studied,  the  tail  region  solver  is  extremely  rapid  and  can  carry 
the  solution  downstream  10  obstacle  nose  radii  in  typically  less  than  10 
CPU  seconds.  As  a  final  run  time  comparison,  we  note  that  the  typical 
CRAY  XMP  CPU  time  required  for  a  single  case  by  the  combined  nose  and 
tail  axisymmetric  model  to  determine  an  analogous  flow  field  is  approxi¬ 
mately  15  seconds. 
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Figure  6  shows  the  calculated  results  for  the  magnetopause  and  bow  wave 

locations  for  conditions  representative  of  those  suggested  by  Slavin  et 

al.  (1985)  as  appropriate  for  Jupiter.  These  conditions  are  that  the 
flow  has  a  free-stream  sonic  Mach  number  Me  of  10  and  ratio  of 

5«o 

specific  heats  y  of  2.0.  and  the  magnetopause  has  an  elliptic  cross 
section  with  constant  ratio  of  major  to  minor  axes  a/b  of  1.75,  corre¬ 
sponding  to  an  eccentricity  e  of  .821.  The  equatorial  magnetopause 
profile  is  taken  to  approximate  that  deduced  from  observations  at 
Jupiter.  Calculation  of  a  converged  solution  for  these  conditions  has 
been  carried  out  over  the  region  extending  from  the  nose  of  the  bow 
shock  to  about  1  magnetopause  nose  distance  or  about  60  Jupiter  radii, 
downstream  of  the  planetary  center.  The  results,  in  the  form  of  traces 
of  the  magnetopause  and  the  bow  shock  in  (a)  each  of  the  ten  equally- 
spaced  azimuthal  grid  planes  from  the  equatorial  plane  to  the  noon- 
midnight  polar  plane  of  symmetry  and  (b)  the  terminator  plane,  display 
that  the  flattening  of  the  bow  shock  is  much  less  than  that  of  the 
magnetopause. 

Properties  of  this  three-dimensional  flow  field  are  presented  in  Figure 
7.  In  that  figure,  contour  plots  of  the  density  p/p^,  pressure  p/p^, 

temperature  T/T  .  flow  speed  V/V  ,  and  the  sonic  Mach  number  Me  in  the 

•O  90  j 

noon-midnight  polar  meridian  and  equatorial  planes  are  presented 
together  with  an  isometric  drawing  of  the  magnetopause  and  bow  shock  in 
the  equatorial  and  noon-midnight  meridian  planes.  These  contours 
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clearly  indicate  the  smoothness  and  highly-converged  nature  of  the 
three-dimensional  solution  for  these  magnetopause  shapes. 

Figure  8  shows  calculated  results  similar  to  those  of  Figure  6  for  MSoo 

12,  y  -  2.0,  and  a/b  -  1.25,  corresponding  to  e  -  0.6,  selected  as 
representative  of  conditions  at  Saturn.  The  equatorial  magnetopause 
profile  is  that  deduced  by  Slavin  et  al .  (1985)  from  the  Pioneer  and 
Voyager  observations  at  Saturn.  As  in  Figure  6,  the  results  show  that 
the  flattening  of  the  trace  of  the  bow  shock  in  the  terminator  plane  is 
much  less  than  that  of  the  magnetopause.  It  is  also  much  less  flattened 
than  the  bow  wave  for  Jupiter,  as  would  b°  expected  in  view  of  the 
substantially  less  flattened  shape  of  the  Saturn  magnetopause. 

Figure  9  presents  a  comparison  of  the  least-squares  fitted  conic  curves 
of  Slavin  et  al .  (1985)  representing  the  Jovian  magnetopause  and  bow 
shock  crossings  observed  near  the  equatorial  plane  with  the  model - 
predicted  bow  shock  shapes  calculated  for  the  same  equatorial  magne¬ 
topause  shape  and  for  four  values  for  a/b.  namely  1.00  (circular  cross 
section).  1.50,  1.75  and  2.00,  that  bracket  the  value  of  5/3  suggested 
by  Slavin  et  al.  (1985).  The  results  show  that  the  bow  shock  is 
farthest  out  for  the  axi symmetric  magnetosphere  and  moves  progressively 
neare-  as  the  flattening  of  the  magnetospheric  cross  section  Increases. 
This  Is  appropriate  because  the  cross-section  area  of  the  magnetosphere 
is  largest  for  the  axisymmetric  case,  and  diminishes  progressively  as 
the  magnetopause  is  flattened  from  a/b  of  1.00  to  2.00,  thereby  present¬ 
ing  less  of  an  obstruction  to  the  flow. 
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Comparison  with  the  observations  for  Jupiter  show  that  they  are  in  close 
accord  with  the  calculated  results  for  a/b  -  1.75,  the  magnetospheric 
shape  of  our  set  of  examples  with  elliptic  cross-section  closest  to  that 
proposed  for  Jupiter  by  Slavin  et  al .  (1985).  This  finding  is  also 
consistent  with  the  result  of  Engle  and  Beard  (1980)  in  which  the  three- 
dimensional  shape  of  the  Jovian  magnetosphere  was  calculated  for  a 
planetary  dipole  field  supplemented  by  an  equatorial  current  sheet 
chosen  to  reproduce  the  observed  Jovian  magnetic  field  near  the  equa¬ 
torial  plane.  All  of  this  confirms  the  general  expectation  that, 
provided  the  solution  is  determined  for  a  suitably  flattened  magne¬ 
topause.  the  gasdynamic  convected  field  model  should  provide  an  even 
better  representation  for  the  magnetosheath  plasma  and  field  properties 
at  Jupiter  than  at  earth  because  of  the  high  sonic  and  Alfvenic  Mach 
numbers  that  prevail  there. 

On  the  other  hand,  the  assumption  that  the  flow  is  steady  and  that  the 
dipole  field  is  stationary  is  less  appropriate  for  Jupiter  and  Saturn 
than  for  the  earth  or  other  terrestrial  planets  because  of  the  immense 
size  of  their  magnetospheres  and  bow  shocks.  Whereas  the  geocentric 

4 

distance  to  the  nose  of  the  earth's  bow  shock  Is  about  8.2x10  km.  the 
corresponding  distances  for  Jupiter  and  Saturn  are  5.3x10®  km  and 
1.75x10®  km.  Note  that  both  of  the  latter  are  greater  than  the  1.4x10® 
km  diameter  of  the  sun. 

Less  than  4  minutes  is  thus  required  for  the  solar  wind  with  a  velocity 
of  400  km/s  to  flow  a  distance  equal  to  that  from  the  nose  of  the  bow 
shock  to  the  terminator  plane  for  the  earth.  For  the  forward  part  of 
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the  magnetopause  and  the  surrounding  flow  field  to  reach  a  quasi -steady 
state,  it  may  be  expected  that  a  period  of  the  order  of.  say.  twice  this 
time,  or  8  minutes.  Is  required.  For  Jupiter  and  Saturn,  the 
corresponding  solar  wind  travel  times  are  about  220  and  75  minutes,  so 
that  the  times  required  to  reach  a  quasi -steady  state  may  be  expected  to 
be  of  the  order  of  7.3  and  2.5  hours,  respectively.  These  are 
exceptionally  long  periods  for  conditions  in  the  solar  wind,  parti¬ 
cularly  the  magnetic  field,  to  remain  essentially  steady. 

A  second  source  of  unsteady  conditions  in  the  Jupiter  magnetosphere  is 
associated  with  the  coning  of  the  magnetic  dipole  axis  as  a  result  of 
planetary  spin  and  the  angular  offset  of  this  axis  from  the  spin  axis. 
These  effects  can  be  disregarded  rather  safely  for  the  earth,  even 
though  there  is  an  angle  of  11.5°  between  these  two  axes,  because  the 
solar  wind  flowing  at  400  km/s  will  travel  1.7xl07  km  2700  earth  radii 
in  12  hours,  a  distance  equal  to  about  200  times  the  distance  from  the 
bow  shock  nose  to  the  terminator  plane.  Furthermore,  it  has  been  long 
established,  see  Spreiter  and  Briggs  (1962)  or  Spreiter  and  Stahara 
(1985)  and  numerous  studies  of  the  observed  positions  of  magnetopause 
crossings,  that  the  shape  of  the  earth's  magnetopause  is  remarkably 
Insensitive  to  the  ±11.5°  daily,  or  even  the  combined  ±35°  yearly- 
daily.  variation  of  the  tilt  of  the  dipole  relative  to  the  direction  of 
the  incident  wind.  The  effects  would  also  appear  to  be  negligible  for 
Saturn  because  Its  dipole  axis  is  aligned  within  1°  of  the  spin  axis  and 
there  is  thus  no  significant  coning  of  the  dipole  axis. 
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The  matter  is  not  so  simply  dismissed  for  Jupiter,  however.  The  dipole 
magnetic  field  for  this  planet,  which  rotates  once  about  every  10  hours. 
Is  inclined  11°  to  the  spin  axis.  Since  a  400  km/s  solar  wind  flows 
about  7.2x10®  km  in  5  hours,  it  follows  that  the  solar  wind  advances 
only  about  1.3  times  the  distance  from  the  bow  shock  nose  to  the 
terminator  plane  during  the  time  the  magnetic  axis  makes  half  a 
revolution  around  the  spin  axis.  This  clearly  would  have  to  be  taken 
into  account  in  the  calculation  of  the  magnetic  field  within  the  magne¬ 
tosphere  and  its  tail,  and  perhaps  the  shape  of  the  magnetopause  down¬ 
stream  of,  say.  the  terminator  plane.  On  the  basis  of  the  insensitivity 
of  the  magnetopause  shape  cited  above  for  the  earth,  we  believe,  how¬ 
ever.  that  the  results  provided  by  the  steady-state  model,  particularly 
those  upstream  of  the  terminator  plane,  are  not  adversely  affected  to 
any  significant  degree  for  many  applications. 

Figure  10  shows  the  corresponding  results  for  conditions  suggested  by 

Slavin  et  al .  (1985)  as  representative  of  Saturn.  These  are  that  the 

flow  has  a  free-stream  sonic  Mach  number  Mc  of  12,  a  value  of  2.0  for 

S<*> 

the  ratio  of  specific  heats  y,  and  is  about  a  magnetopause  having  the 
indicated  shape  in  the  equatorial  plane  and  an  elliptic  cross  section 
with  constant  ratio  of  major  to  minor  axes  a/b  somewhere  between  that 
for  the  earth  and  Jupiter.  Calculated  equatorial  plane  traces  of  the 
bow  shock  are  presented  for  magnetopause  shapes  with  a/b  -  1.00 
(circular  cross  section),  1.15,  1.25,  1.35  and  1.45  together  with  the 
least-square  curve  fitted  to  the  observations  at  Saturn  by  Slavin  et  ai. 
(1985),  and  also  the  new  least-square  conic  curve  included  on  Figure  2. 
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The  results  calculated  for  a/b  -  1.25  provide  a  good  match  to  the  least- 
square  fitted  conic  curves  near  the  nose  of  the  bow  shock.  However,  the 
latter  flares  out  much  farther  from  the  planet  than  not  only  the  calcu¬ 
lated  results  for  a/b  -  1.25,  but  the  results  for  all  the  calculated 
cases,  including  the  axisymmetric  magnetopause.  This  result  is  surpris¬ 
ing  in  view  of  the  generally  satisfactory  agreement  between  the  results 
of  the  gasdynamic  calculations  and  the  observations  for  all  the  other 
planets,  and  has  been  commented  on  by  Russell  (1985)  and  Slavin  et  al . 
(1985). 

Prompted  by  suggestions  of  Dr.  R.  P.  Lepping  of  NASA  Goddard  Space 
Flight  Center,  we  have  re-examined  the  procedures  used  by  Slavin  et  al . 
(1985),  to  determine  their  least-square  conic  curves  for  the  bow  shock 
and  magnetopause  for  Saturn  and  arrived  at  the  new  least-square  conic 
curve  included  in  Figures  2  and  10.  This  curve  is  fitted  to  the  stagna¬ 
tion-pressure-corrected  observations  using  the  same  techniques  as 
employed  by  Slavin  et  al .  (1985),  but  including  the  data  from  Voyager  1. 
Those  observations  were  not  included  in  the  original  determination 
because  the  shock  crossings  occurred  approximately  two  magnetopause  nose 
distances  downstream  of  the  terminator  plane,  and  did  not  meet  the 
general  modeling  criteria  set  up  and  used  previously  by  Slavin  and 
Holzer  (1981)  and  Slavin  et  al .  (1983)  In  analyses  of  observations  of 
the  terrestrial  planets.  To  ensure  comparability  of  the  models  for  the 
magnetopause  and  bow  shock  shapes  for  Jupiter  and  Saturn  with  similar 
ones  for  Venus,  earth  and  Mars,  the  Voyager  1  crossings  were  therefore 
not  used  in  determining  the  least-square  conic  curves  of  Slavin  et  al. 
(1985).  Since  comparisons  of  observations  of  the  terrestrial  and  Jovian 
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planets  is  not  an  objective  for  this  study,  such  a  reason  is  not  so 
compelling  here.  Moreover,  the  location  of  the  Voyager  1  outbound 
crossing  shown  in  Figure  1  suggests  that  its  inclusion  might  signifi¬ 
cantly  reduce  the  large  flare  of  the  least-square  curve  to  the  Saturn 
shock  and  improve  the  agreement  with  both  the  calculated  results  and 
also  the  shock  shapes  observed  at  the  other  planets.  The  results 
presented  in  Figure  10  show  that  this  expectation  is  realized  to  only  a 
limited  degree,  however,  since  even  the  new  least-square  conic  curve 
flares  out  beyond  the  calculated  curves  by  the  terminator  plane. 

Both  the  original  and  the  new  least-square  curves  are  conic  sections 

described  in  terms  of  nondimensional  cylindrical  coordinates  X/R  ,  R/R 

oo 

by  equations  of  the  form  X  -  XQ  +  rcos0  and  R  -  rsine,  where  r  -  L/(l  + 
ccose  ).  RQ  -  19.7  x  (2.1  x  10'10/3.2  x  10'10)1/6  -  18.36,  and  0  is  the 

angle  measured  from  the  x  axis  with  0  -  0  at  the  magnetosphere  nose. 

The  original  curve  determined  without  use  of  the  outbound  Voyager  1 
observations  is  defined  by  L  -  55.4,  XQ  -  6,  and  e  -  1.71.  The 

corresponding  values  for  the  new  curve  are  L  -  51.8,  XQ  -  5.  and  c  - 

1.39;  and  those  for  the  magnetopause  are  L  -  30.8,  XQ  -  5.  e  -  1.09,  and 

Rq  -  19.7.  In  neither  case  was  any  attempt  made  to  eliminate  Voyager  2 

data  on  the  basis  of  it  being  atypical  because  of  Saturn  possibly  being 
in  the  magnetic  tail  of  Jupiter,  although  It  is  evident  from  Figure  1 
that  elimination  of  those  data  would  lead  to  a  least-square  fitted  conic 
curve  that  would  not  only  be  in  improved  agreement  with  our  calculated 
results  but  with  general  experience  for  all  the  other  planets  for  which 
such  studies  have  been  made.  It  should  .be  noted  that  the  procedure 
employed  in  determining  these  least-square  curves  is  increasingly 
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suspect  with  increasing  distance  from  the  nose  because  no  consideration 
is  given  in  their  determination  to  matching  the  asymptotic  directions  of 
the  bow  shock.  This  would  not  pose  a  problem  if  adequate  data  were 
available  for  downstream  shock  crossings,  but  that  Is  not  the  situation 
for  Saturn. 

Although  the  stagnation  pressure  correction  introduced  through  Rq 

accounts  for  a  significant  part  of  the  inward  location  of  the  Voyager  1 
outbound  shock  location,  the  overall  result  is  a  significant  reduction 
in  the  amount  of  flaring  at  the  flanks  of  the  bow  shock.  Even  so,  down¬ 
stream  of  the  terminator  plane,  the  revised  least-square  conic  curve 
extends  outside  the  location  of  the  calculated  shock  for  an  axisymmetric 
magnetosphere  having  the  same  shape  as  the  least-square  conic  curve 
representing  the  low-latitude  magnetopause  of  Saturn.  It  would  appear 
from  these  comparisons  that  magnetospheric  flattening  for  Saturn  may  be 
about  one-third  that  for  Jupiter,  but  the  differences  in  the  shock  loca¬ 
tions  near  the  terminator  plane  require  further  examination. 

Concluding  Remarks 

The  development  of  a  three-dimensional  computational  model  for 
determining  the  solution  of  the  gasdynamic  portion  of  the  gasdynamic 
convected  magnetic  field  model  for  solar  wind  flow  past  a  magne- 
tolonopause  of  nonaxlsymmetric  shape  is  described.  Specific  application 
of  the  model  is  made  to  Jupiter  and  Saturn  where,  as  inferred  from 
observations,  significant  equatorial  broadening  of  the  magnetopause 
cross  sections  exists.  Results  from  the  new  three-dimensional  model 
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converge  continuously  to  those  calculated  previously  for  axisymmetric 
magnetoionospheres  that  were  determined  using  completely  different 
computational  procedures,  thereby  verifying  the  accuracy  of  both 
methods. 

Because  the  cost  of  determining  these  three-dimensional  solutions  is  two 
orders  of  magnitude  greater  than  for  axisymmetric  flows,  and  the 
differences  between  the  results  are  small  for  small  departures  from 
axisymmetry.  we  suggest  that  the  axisymmetric  model  be  used  whenever  the 
magnetoionopause  cross  section  is  nearly  circular,  and  that  the  new  non- 
axisymmetric  model  be  employed  only  for  those  cases  in  which  substantial 
departures  from  axisymmetry  exist.  Examples  of  the  latter  are  Jupiter 
and  Saturn,  for  which  previous  comparisons  between  the  observations  and 
results  calculated  for  axisymmetric  magnetospheres  were  shown  to  be  in 
substantially  poorer  agreement  with  observations  than  for  the  terres¬ 
trial  planets. 

Comparisons  of  results  from  the  new  model  for  a  family  of  three- 
dimensional  nonaxisymmetric  magnetospheres  having  elliptic  cross 
sections  of  various  eccentricities  with  the  observed  magnetopause  and 
bow  shock  crossings  at  Jupiter  and  Saturn,  all  of  which  occurred  at  low 
latitudes,  provide  a  new  means  to  estimate  the  amount  of  equatorial 
broadening  of  the  magnetospheres  of  these  planets.  It  is  found  that 
essentially  perfect  agreement  with  the  observations  is  obtained  for 
Jupiter  when  a  value  of  1.75  is  used  for  the  ratio  a/b  of  major  to  minor 
axes  of  the  cross  section.  For  Saturn,  similarly  good  agreement  with 
the  observations  is  found  in  the  nose  region  when  a  value  of  1.25  is 
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used  for  a/b.  Both  of  these  values  are  in  essential  agreement  with  pre¬ 
vious  estimates  of  5/3  and  4/3  for  a/b  based  on  other  considerations. 

Near  and  downstream  of  the  terminator  of  Saturn,  however,  the  least- 
square  fitted  conic  curve  (Slavin  et  al.,  1985)  used  to  represent  the 
shape  of  the  bow  shock  flares  out  considerably  more  than  Indicated  by 
any  of  the  model  predictions.  Including  that  for  an  axisymmetric  magne¬ 
topause.  In  an  attempt  to  gain  a  better  understanding  of  the  reason  for 
this  behavior,  we  have  reexamined  the  fitting  procedure  applied  to  the 
Saturn  data  and  produced  a  new  least-square  conic  curve  that  includes 
the  observations  of  Voyager  1  spacecraft  not  utilized  in  the  original 
determination.  Those  observations  were  not  used  in  the  original 
determination  because  the  bow  shock  crossings  were  too  far  downstream  of 
the  planet  to  meet  criteria  set  previously  in  studies  of  the  terrestrial 
planets.  The  new  curve  matches  the  calculated  shape  of  the  bow  shock 
for  a/b  -  1.25  for  a  greater  distance  from  the  nose  of  the  bow  shock 
than  the  original  curve,  but  also  flares  out  more  than  the  cal  ulated 
results  at  and  downstream  of  the  terminator.  An  inescapable  conclusion 
of  these  comparisons  of  calculated  and  observed  shock  shapes  is  that  the 
currently  available  observational  data  for  bow  shock  crossings  at  Saturn 
are  too  few  in  number  and  too  scattered  to  provide  a  well-defined 
representation  of  the  shape  of  the  bow  shock  in  the  terminator  region. 
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FIGURE  Till-S 


Figure  1.  Locations  of  the  Jupiter  and  Saturn  bow  shock  and 
magnetopause  crossings  observed  by  Pioneer  10.  11  and  Voyager  1.  2 
spacecraft  (Slavln  et  al..  1985). 

Figure  2.  Least-square  fitted  conic  curves  representing  the 
observations  of  the  Jovian  and  Saturnian  bow  shock  and  magnetopause  in 
the  planetary  equatorial  plane  together  with  predictions  (dashed  lines) 
of  the  gasdynamic  bow  shock  model  for  an  axisymmetric  magnetopause 
matching  the  fitted  magnetopause  curve  near  the  equatorial  plane. 

Figure  3.  A  conceptual  representation  of  magnetopause  cross  sectional 
shapes  based  upon  observations,  theoretical  magnetic  field  models,  and 
gasdynamic  modeling  results. 

Figure  4.  Illustration  of  flow  field  subdivision  and  coupled  flow 
solver  combination  incorporated  in  the  three-dimensional  gasdynamic 
solar  wind  planetary  interaction  model.. 

Figure  5.  Perspective  cutout  view  of  one  quadrant  of  a  typical 
computational  grid  employed  In  the  three-dimensional  model. 

Figure  6.  Predicted  bow  shock  and  corresponding  magnetopause  locations 
In  the  terminator  plane  and  in  the  10  equally-spaced  azimuthal  planes 
employed  in  the  three-dimensional  model  for  an  equatorial ly  broadened 
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Jovian  magnetopause  having  an  elliptic  cross  section  with  a/b  -  1.75  for 
Mc  -  10  and  y  -  2.0. 

o®° 

Figure  7.  Contours  of  density,  pressure,  temperature,  velocity 

magnitude  normalized  by  free  stream  values,  and  Hach  number  In  the 

equatorial  and  polar  meridian  planes  for  the  Jovian  magnetopause  with 
a/b  -  1.75,  Mc  -  10  and  y  -  2.0. 

(a)  Noon-midnight  polar  meridian  plane 

(b)  Equatorial  plane 

Figure  8.  Predicted  bow  shock  and  corresponding  magnetopause  locations 

in  the  terminator  plane  and  in  the  10  equally-spaced  azimuthal  planes 

employed  in  the  three-dimensional  model  for  an  equatorially  broadened 

Saturnian  magnetopause  having  an  elliptic  cross  section  with  a/b  -  1.25 
for  Mc  -  12  and  y  -  2.0. 
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Figure  9.  Comparison  of  a  least-square  fitted  conic  curve  representing 
the  observational  data  for  the  Jovian  bow  shock  In  the  planetary 
equatorial  plane  with  three-dimensional  bow  shock  predictions  for  magne¬ 
topause  shapes  having  elliptic  cross  sections  with  different  ratios  a/b 

of  major  to  minor  axes  representing  various  amounts  of  equatorial 
broadening  for  Ms#o  -  10  and  y  -  2.0. 

Figure  10.  Comparison  of  a  least-square  fitted  conic  curve  representing 

the  observational  data  for  the  Saturnian  bow  shock  in  the  planetary 

equatorial  plane  with  three-dimensional  bow  shock  predictions  for 

magnetopause  shapes  having  elliptic  cross  sections  with  different  ratios 

a/b  of  major  to  minor  axes  representing  various  amounts  of  equatorial 
broadening  for  Mc  -  12  and  y  -  2.0. 
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Figure  1.  Locations  of  the  Jupiter  and  Saturn  bow  shock  and 
•agnetopause  crossings  observed  by  Pioneer  10,  11 
and  Voyager  1,  2  spacecraft  (Slavln  et  al.,  1965). 
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Figure  2. 


Least-square  fitted  conic  curves  representing  the 
observations  of  the  Jovian  and  Saturnian  bow  shock 
and  Magnetopause  In  the  planetary  equatorial  plane 
together  with  predictions  (dashed  lines)  of  the 
gasdynamic  bow  shock  Model  for  an  axisyMmetric  Mag¬ 
netopause  Matching  the  fitted  Magnetopause  curve 
near  the  equatorial  plane. 
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Flqure  3.  A  conceptual  representation  of  magnetopause  cross 
sectional  shapes  based  upon  observations,  theoreti¬ 
cal  magnetic  field  models,  and  gasdynamic  modeling 
results.  * 
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Figure  4.  Illustration  of  fl ow  field  subdivision  and  coupled 
flow  solver  combination  incorporated  In  the  three- 
dimensional  gasdynamic  solar  wind  planetary  Inter¬ 
action  model. 
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Figure  5.  Perspective  cutout  view  of  one  quadrant  of 
typical  computational  grid  employed  In  the  three 
dimensional  model. 
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Figure  6.  Predicted  bow  shock  and  corresponding  magnetopause 
locations  in  the  terminator  plane  and  in  the  10 
equally-spaced  azimuthal  planes  employed  in  the 
three-dimensional  model  for  an  equatorial ly  broad¬ 
ened  Jovian  magnetopause  having  an  elliptic  cross 
section  with  a/b  ■  1.75  for  MSoc  ■  10  and  y  »  2.0. 
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(b)  Equatorial  plane 


Figure  7.  Contours  of  density,  pressure,  temperature,  velo 
city  magnitude  normalized  by  free  stream  values 
•nd  Mach  number  In  the  equatorial  and  polar  meri 
iian  planes  for  tne  Jovian  magnetopause  with  a/b 
*•75,  •  10  and  y  •  2.0. 
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Figure  9.  Comparison  of  a  least-square  fitted  conic  curve 
representing  the  observational  data  for  the  Jovian 
bow  shock  in  the  planetary  equatorial  plane  with 
three-dimensional  bow  shock  predictions  for  magne¬ 
topause  shapes  having  elliptic  cross  sections  with 
different  ratios  a/b  of  major  to  minor  axes  repre¬ 
senting  various  amounts  of  equatorial  broadening 
for  MSoo  ■  10  and  y  ■  2.0. 
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Figure  10.  Comparison  of  a  least-square  fitted  conic  curve 
representing  the  observational  data  for  the  Satur¬ 
nian  bow  shock  In  the  planetary  equatorial  plane 
with  three-dimensional  bow  shock  predictions  for 
magnetopause  shapes  having  elliptic  cross  sections 
with  different  ratios  a/b  of  major  to  minor  axes 
representing  various  amounts  of  equatorial  broad¬ 
ening  for  HSflo  ■  12  and  y  ■  2.0. 
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